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1. Introduction 


Computation of the flow field inside a space shuttle main engine (SSME) requires 
the application of the state-of-the-art CFD technology. Several computer codes are under 
development to solve three dimensional Navier— Stokes equations for analyzing the SSME 
internal flow, such as the flow through the hot gas manifold. The computational methods 
used in the Navier— Stokes codes fall into two major categories: finite difference and finite 
element methods. Some of the algorithms are designed to solve the unsteady compress- 
ible Navier— Stokes equations, either by explicit or by implicit factorization methods, using 
several hundred or thousands of time steps to reach a steady-state solution asymptotically. 
Other algorithms attempt to solve the steady-state equations by relaxation methods. All 
of them require body-fitting curvilinear grids with sufficient resolution. Grid requirements, 
however, differ greatly with the region being modelled and the algorithm used. Implicit 
factorization based on finite differences typically uses global numerical transformations 
whereby the transformed grid in the computational space is uniform and rectilinear. This 
requires the grid to have indices which are separable in the three directions for three di- 
mensional problems, and also be reasonably smooth. However, such requirements may 
introduce grid singularities when complicated domains are discretized. Flow solver algo- 
rithm will have to deal with such grid singularities. Explicit schemes and finite element al- 
gorithms have less stringent requirements on the grid structure. However, explicit schemes 
are slow to converge because of the stability limitations on time step, particularly for large 
scale viscous problems. 

The finite element method is characterized by three basic features which are credited 
for the enormous success the method has enjoyed in the solution of practical engineering 
problems. The first feature is that every computational domain is viewed as a collection of 
simple subdomains, called finite elements. This feature allows us to represent complicated 
geometries as assemblages of simple parts. It is a desirable feature in the solution of flow 
problems in complex configurations, not only to describe the complex geometry but also to 
choose the most suitable computational grid for a particular flow. This feature also allows 
us to place or remove any obstructions routinely into the flow field. The second feature is 
that over each element the solution is represented by polynomials of desired degree. This 
allows us to compute the solution as a continuous function of position instead of at selected 
few points. The third feature is that the relationship (i.e., the algebraic equations) between 
the solution and its dual variables is developed using a variational method, such as the 
Galerkin method. The boundary conditions are then applied on the algebraic equations 
directly before solving. The three features of the finite element method also allow the easy 
development and interfacing of pre- and post-processors, and user-defined subroutines for 
equations for state and turbulence models. 

The Galerkin finite element method (i.e., the weight functions are the same as the 
approximation functions) applied to flow problems always results in implicit schemes. The 
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weighted-residual (or Petrov-Galerkin) method, in which the weight functions are differ- 
ent from the approximation functions, can be used in conjunction with explicit schemes 
to obtain explicit final equations. For example, by selecting the weight functions to be 
orthogonal to the approximation functions, the mass matrix can be diagonalized. How- 
ever, such considerations are entirely in the interest of obtaining explicit schemes and not 
necessarily in the interest of accuracy or even computational efficiency. In the current 
project an implicit finite element scheme with suitable dissipation terms for stability is 
developed. A relaxation procedure, known as the locally implicit scheme is developed to 
solve the coupled set of algebraic equations efficiently. 

Allowing the possibility of unstructured grids is important for discretizing complex 
flow domains efficiently and also for adding the features of solution-adaptive grids. For 
grids with large numbers of nodes, direct solution procedures for the finite element equa- 
tions become impractical. Thus we have undertaken the development of a new iterative 
algorithm for the solution of implicit finite element equations without assembling global 
matrices. It is an efficient iteration scheme based on a modified non-linear Gauss-Seidel 
iteration with symmetric sweeps. This algorithm is analyzed for a model equation and is 
shown to be unconditionally stable. This analysis is reported in the next Section. 

The locally implicit scheme is unconditionally stable based on local linearized anal- 
ysis. However, for strongly convective flows there is a possibility of non-linear numerical 
instabilities occurring in some parts of the flow domain and eventually destabilizing the 
entire flow domain. We have added adaptive artificial dissipation terms of third order to 
the finite element approximations similar to Jameson and others^ 1 ). These are designed 
to suppress non-linear instabilities if they appear and at the same time be much smaller 
than the real viscosity terms in viscous zones. 

In numerical schemes for solving fluid flow equations, there is some degree of uncer- 
tainty as to the imposition of boundary conditions on some of the variables at different 
types of boundaries, particularly at the inflow and outflow boundaries. In the current fi- 
nite element code we have developed special procedures to compute the required flux terms 
at the boundary surfaces to the same degrees of accuracy as in the interior. We expect 
that our technique of computing the required surface fluxes iteratively, together with the 
interior flow variables, should minimi ze the uncertainties in the imposition of boundary 
conditions. 

The locally implicit scheme is tested on a variety of problems. It has been shown 
to be efficient with multi-grid acceleration procedures for elliptic problems by Reddy and 
Nayani* 2 ) and for inviscid compressible flows from transonic to supersonic Mach numbers 
by Reddy and Jacocks (3 >. Reddy, Reddy and Nayani (4) have developed this scheme for 
viscous flow problems. We developed a 2-D test code for solving unsteady compressible 
Navier-Stokes equations with finite volume approximation, which is a special case of the 
finite element approximation. This code has been used to check various features of the 
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locally implicit solution algorithm. We have also added an algebraic turbulence model 
developed by Baldwin and Lomax^ 5 ). 

Results for a series of test problems are presented in this report. The finite element 
code has been tested for Couette flow, described in Schlichting^ , which is a flow under a 
pressure gradient between two parallel plates in relative motion. Another problem that has 
been solved is viscous laminar flow over a flat plate. As a test case for the locally implicit 
scheme, the 2-D finite volume code has been applied to compute subsonic and transonic 
viscous flows over airfoils for both laminar and turbulent cases. The general 3-D finite 
element code has been used to compute the flow in an axisymmetric turnaround duct at 
low Mach numbers. 
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2. Locally Implicit Scheme for a Model Equation 

Locally implicit scheme is a relaxation method for solving the non-linear finite element 
equations approximating the Navier-Stokes equations. It is a point iteration method at 
each time step. However, it is not necessary for the iteration to converge fully at each time 
step if we are interested in computing the time asymptotic steady-state solutions. The 
analysis of the consistency, stability and hence convergence of the scheme is presented for 
a model equation for the Navier-Stokes equations. 

Consider a one-dimensional convection-diffusion equation, 

du du d 2 u 

dt °dx V dx 2 

Finite element approximation at a node j on a unif orm mesh for equation (2.1) can be 
written as 

Ft J u *’ iz + / + "t:) lt dx = 0 < 2 - 2 ) 

where <f>j is a global test function corresponding to the node j. For a linear element 
approximation, equation (2.2) gives 


( 2 . 1 ) 


f {ib- + b + b +i } + (2b) (u ' +i > 

~ (^ 2 ) “ 2u J + u i+i) = 0 


Implicit time integration gives 


J Au '- 1 + | Au > + + f W? - u i-i) 

- R - 2 “" + ’ + “"« ) = 0 


(2.3) 


(2.4) 


where Auj = tt" +1 — it" 


C = aAt/Ax , R = uAt/Ax 2 

Equation (2.4), together with appropriate boundary conditions, gives a system of linear 
equations which can be solved easily in one-dimension and this scheme is unconditionally 
stable. However, the system of equations becomes too large in multi-dimensions and vari- 
ous types of sparse matrix solvers are developed in the literature, but they are usable only 
with a modest number of nodes. Alternately, we develop a relaxation scheme to solve (2.4) 
approximately at each time step. The scheme is a modification of the symmetric Gauss- 
Seidel iteration. The basic Gauss— Seidel iteration, even with symmetric sweeps, is unstable 
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for a whole range of Courant number C in equation (2.4). The present modification makes 
it unconditionally stable. Rewrite the equation (2.4) in delta form as 

+ g A “'+> + f (Axy+l - (2.5) 

— R(Auj-\ — 2A tij + Auj+j) = Res" 

where 

- -f W+i - <*"-.) + R (“"-i - + “’+i) (2 6) 

As Auj = u" +1 — u” — ► 0 as n — + oo, we obtain the asymptotic steady-state solution as the 
Resj function is driven to zero. This process may be speeded up and made more robust 
by choosing a value for R on the left side of equation (2.5) larger than the value of R on 
the right side of equation (2.5). To analyze this process we use the notation R for R on 
the left side of equation (2.5). It may be noted that we can always obtain time accurate 
solution, if that is required, by choosing R = R. We solve for Auj at each time step by a 
modified Gauss-Seidel iteration: 

Au$ m+1) = A u$ m) + duj, A u$ 0) = 0 (2.7) 


Left-to-right sweep yields 


-duj + -du j+1 + — du j+1 
— R(—2duj + duj+i) = RES 


( 2 . 8 ) 


where 


RES =R C >^ - 

+f (A-#! - Au£r>) - R (Au^" - 2Au(”> + 

Now we approximate duj+i a duj and replace C by its absolute value \C\ on the left side 
of equation (2.8), to accommodate convection velocity direction either in or opposite to 
the iteration sweep direction. This leads to an explicit expression for duj: 

+ du, = RBS (2.10) 

Right-to-left sweep is defined similarly. A symmetric iteration sweep consists of a left-to- 
right sweep followed by a right-to-left sweep. It may be noted that duj is the iterative 
correction to the time change iterates Au^ and if the iteration process is convergent, 
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RES — ► 0 and the equation (2.5) can be satisfied as accurately as we wish by carrying 
out the necessary number of symmetric iteration sweeps. The approximations made in 
the iteration do not affect the actual solution itself. Thus the iteration equations are 
consistent with the basic equations. One or two symmetric sweeps per time step are usually 
sufficient for obtaining steady-state solutions. Local stability analysis can be carried out 
by computing the amplification factor of discrete Fourier modal solutions per time step. 
In this analysis, we seek modal solutions of the equations (2.9) and (2.10) in the form 

u? = v n e iJ( , 0 < £ = aAx < ic 
Au$ m) = At> (m V*, m = 0 , 1 , . . . 

u " +1 = v n+1 e ij * 

For a single symmetric sweep per time step (m = 0, 1), 

n n+1 = v n + Av (2) = g(t)v n 


where </(£) is known as the amplification factor from one time step to the next and is given 
by 

9(0 = 1 + T- 1 + T - , 0 < £ < ? r 

/l3 L «1 . 

r = — Ci sin £ +2J?(cos £ — 1) 


+ f+S-i) 


b 


5 

6 



( 2 . 11 ) 


A necessary condition for stability is |<j(£)| ^ 1- It can be shown that |p(£)l * 8 indeed < 1 
unconditionally. It is also desirable to have |^(£)| < 1 as much as possible for £ closer to n 
which represents the range of high frequency modes of the solution. Figure 1 shows plots of 
|y(£)| versus £ for different Courant numbers for R = R = Figure 2 shows plots of |g| 
versus £ for C = 10, R = R and R takes different values. Figure 3 shows the plots for C = 
10, R = 2 R and R takes different values. Numerical plots of |y| against £ confirm that the 
scheme is unconditionally stable. However, very large Courant numbers are not necessarily 
the best. Courant number C ~ 10 and R = 2 R — ► 4 R seem desirable ranges. Amplification 
factors corresponding to two or more symmetric modified Gauss-Seidel iterations have 
similar behavior. Thus we establish unconditional stability for the modified Gauss-Seidel 
iteration scheme for the convection-diffusion equation. Similar stability can be shown 
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when the diffusion term is replaced by a 4th difference term of the type that is used as 
artificial viscosity term of third order for suppressing non-linear instabilities for convection 
dominated flows. It is possible to use artificial viscosity terms which are smaller than the 
truncation terms of the second order accurate finite element approximations. In the present 
Navier-Stokes finite element code where we compute all terms to full second order accuracy, 
artificial dissipation terms, which are an order of magnitude smaller then truncation error, 
are included to suppress non-linear instabilities. Stability analysis of the model equation 
indicates that the locally implicit scheme is unconditionally stable in a local linearized 
sense. 
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3. Locally Implicit Scheme for Navier— Stokes Equations 


Many algorithms designed to solve the unsteady compressible Navier-Stokes equations 
use either explicit methods or implicit factorization methods. Finite element approxima- 
tions usually yield implicit equations. These are solved by explicit time integration meth- 
ods after making additional approximations. Explicit methods may take thousands of 
time steps to converge. Solving them implicitly with Newton iteration is possible, but the 
matrix storage requirements for the resulting algebraic equations and the solution process 
make it prohibitive even for modest size three dimensional flow problems. There are other 
algorithms based on relaxation methods. We have developed a locally implicit method for 
solving the non-linear finite element approximations for 3-D Navier-Stokes equations at 
each time step. 

The method is based on a relaxation procedure for solving the finite element equations 
corresponding to each node iteratively. The equations for the elements surrounding a 
particular node are evaluated based on the latest iterates for the flow variables at the 
nodes around it and the solution is updated at that node by a modified Gauss-Seidel 
iteration. This procedure does not require the assembly of a global matrix, in contrast 
to the standard finite element algorithms. It does not require the solution of a system 
of large number of equations. Thus it is a matrix-free implicit finite element algorithm. 
An additional feature of the algorithm is that while it uses tri-linear approximations for 
the flow variables in quadilateral (brick) elements, all the non-linear fluxes in the Navier- 
Stokes equations are evaluated without any further linear approximation. The fluxes are 
non-linear and are computed accordingly. This assures the second order spatial accuracy 
of the scheme even for unstructured grids. 

3.1 Finite Element Approximations 

The unsteady, compressible Navier-Stokes equations are written in conservation form 
as 

{^}+^{f'}+V(F'} = {°} (3.1) 

where 



{F 7 } and {F v } represent the inviscid and viscous fluxes respectively. Details of these 
equations are given in Appendix I. 
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The variational form (weak form) of equation (3.1) over an element f l e is written as 


0 = L ( { * }T {fr} ~ ' {p ’ + ^“>) dV + jf {*) T {F„}d5 (3.2) 


where {$} are test functions. They are tri-linear functions for linear finite element ap- 
proximation and piecewise constants for finite volume approximations. F n = (F v + F 1 ) • n 
where n is the outward drawn unit normal to the surface S e of the element f 1*. The con- 
servation variables U = (U Q ,a = 1, • • • 5) are approximated by the interpolation functions 
’i ' j as 
N 

K.-E 0'«/(*,v,*)b{»}{K.} (3.3) 

>=1 

where 

{*} = {*,*2 •••*„}, {&„} = (cj,t£, ■■■ o?) T 


U£ is the numerical value of the ath component of the conservation flow variable U at 
jth local node of the element fl e . The interpolation functions and test functions $ are 
chosen to be the same for compressible flow equations. N = 8 for tri-linear approximations 
on quadrilateral brick elements. These approximations are done according to the standard 
finite element approximations (Ref. 7). 


Define the total nodal vector of the 


conservation variables at the nodes of an element as 



{U,}) 


rW 

{0} 

{0} 

{0} 

{0}1 


{U 2 } 


{0} 

m 

{0} 

{0} 

{0} 

{to = • 

K * J 

• 

. ; (*]• = 

{0} 

{0} 

{*} 

{0} 

{0} 

5JV x 1 


5 x 5N 

{0} 

{0} 

{0} 

w 

{0} 


l {^5} , 


L {0} 

{0} 

{0} 

{0} 

mJ 


Then 

{tO = I V * | - [*]*{{?}• 

5x1 U.J 

Now the variational statement (2) can be written as 

(0) = jf _ ([*] T |*|{&} - l^*] 7 ' • (/)) dV + / [#] t {F„}<JS 


(3.4) 


(3.5) 


It should be noted at this point that $ and F n are non-linear functions of U and thus 
the integrals involving them can be expressed analytically in terms of the components 
of U . These expressions are long but they can be programmed into the computer code 
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efficiently. The coupled non-linear differential equations (3.5) are discretized in time by 
the Euler implicit scheme as follows: 

±[M e ]{AU e } + {77*}"* +1 = {0} (3.6) 

where 

AU e = (U e ) m+1 - (t/‘) m , m - time level 

[Af *] = / [*] r [tf]dF (3.7) 

Jn « 

{77*} = - f [V*] T -{F}dV+ l [9] T {F n }dS (3.8) 

Jo* Js • 

Details of the expression {77*} in equation (3.8) are given in Appendix II. In the standard 
finite element algorithms, the element equations (3.6) are linearized, usually by Newton 
method, and all the element equations are assembled to derive a global system of linear 
equations which are solved by sparse matrix solvers. For large scale problems the matrices 
involved become too big to be practical. Here we develop a matrix-free relaxation method 
to solve the non-linear equations directly by a modified Gauss-Seidel iteration. 


3.2 Locally Implicit Scheme 


We wish to solve the non-linear finite element equations iteratively at a node i. We 
assume the nodal values of the solution at all the surrounding nodes from their latest 
iterates. The test function 'J'i, corresponding to the node t, in equation (3.6) gives the 
contribution of element fl* to the node i in the finite element approximation. Adding 
similar equations from all the elements surrounding a node ND yields the nodal finite 
element equation. Thus the equations corresponding to a single node, ND are 

£ Al/<} + (KT +1 ) kd = 0 (3.9) 


where U e is replaced by U e for convenience. Thus U e is the conservation variable vector 
at all the nodes of the element e, and the summation in equation (3.9) is over all elements 
e surrounding the node ND. Equation (3.9) represents 5 equations at ND corresponding 
to each of the 5 conservation equations. The ath conservation equation at ND can be 
written as 



<f 

Jd O' 


e 

(ND)" 


ndS 


= 0 


(3.10) 
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where ^(/vd) = ^th * corresponding to the local index of the global node ND in element 
e. For all interior nodes ND, the surface flux integral in equation (3.10) vanishes. This 
equation couples U at all the nodes surrounding the node ND. We develop a modified 
symmetric non-linear Gauss-Seidel iteration to solve the coupled system of non- linear 
equations directly without linearization. This leads to a matrix-free algorithm for the 
solution. 

For a particular time step n, the iteration is carried out as follows. During the iteration 
process, we assume that all t/’s in the ath equation other than U a are known from the 
previous step of the iteration. We solve for AU a at node ND approximately by a modified 
Gauss-Seidel iteration. 

A£/<” +1) = + dU a ,j (3.11) 

for all nodes j where (m + l)th iterates are not available. 

F a(n+1) ~ F a (U n + (3.12) 


at nodes where At/^ m+1 ) is available. At other nodes where only A is available, 

jtor(n+l) ^ p a fan + A tf("0 + 


- F° (U n + A + ^jdU 


(3.13) 


dF 


gplnvis gfVis 


The Jacobian matrices have inviscid and viscous parts — — respectively. 

The inviscid part is approximated by the spectral radii of the Jacobian matrices multiplied 
by identity matrices. 


gplnvis 

— OU ♦(|tz| + a, |v| + a, H + a)/ = 5)2 (3.14) 

where u,v,w are velocity components and a is the speed of sound. The viscous parts of 
the Jacobian matrices are not altered. For the iterative corrections dU's we make the 
approximation, 

dU a j ~ dU Q> (ND) (3.15) 

for all the nodes j at which the latest iterates are not available. dU a< (SD) = dU a ^ where 
* is the local index corresponding to the global node ND. With this approximation, we 
obtain explicit scalar expression for the iterative correction at the node ND, dU a< (SD). 

C dU ai ND = — (3.16) 
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where 



- / ' F‘"'W + / *Ind)F“ M ■ ndS 

J n« Jan « 

The superscript (*) corresponds to the iteration level (m) or (m + 1) which ever is available 
at the nodes surrounding the node ( ND ). 


C = 
+ 



V‘*l ND) IND(j)dV j 
I ■ SR , D) dV + Y, 



■'Z.imj) 


gpaVis 


dV 


(3.18) 


{ 1 , for nodes j at iteration level m 

(3.19) 

0 , for nodes j at iteration level m + 1 

The absolute value sign | • | in the middle integral indicates the absolute values of each of its 
components. In defining the coefficient C, contributions of surface integrals do not exist for 
all interior nodes and they are ignored for boundary nodes for simplicity. Approximations 
made in C to simplify the algorithm while preserving numerical stability for large Courant 
numbers, do not affect the solution which is obtained by driving Res function to zero. One 
iteration sweep starting from the initial node to the final node followed by a reverse sweep 
makes one symmetric sweep. Typically two symmetric sweeps per time step are sufficient 
for obtaining time asymptotic solutions. 


3.3 Surface Flux Computation 

Volume integrals over quadrilateral brick elements are computed by isoparametric 
transformations to a standard cube and by the use of two point Gaussian integration in 
each direction. The details of such computations are available in many books on finite 
element methods. Surface flux computation, however, is less known and the basic idea is 
outlined below. 

Suppose £, 77, ( are the local coordinates and z, y, 2 are global coordinates and we wish 
to compute the surface flux on the surface ( = 1 of an element. 
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Q (x(£ + Ae.u.O. 



*(£,»? + Arj.f)) 

<f F ■ ndS = <f FdS (3.20) 

A=i •/<=! 

d5 = nrfS = 0>x0g 

= (x c A£, y^A£, z* A£) x (x,A» 7 , y. A t/, x, At;) (3.21) 

„ f ^(y^) d(z,x) ^(x,y) \ 

U^)* ace,*? • 

® F • dS can now be computed with Gaussian integration in ( and rj directions, at C 
J(=i 

= 1. The values of F and the surface Jacobians are evaluated at the Gaussian points 
on the surfaces of the elements, in contrast to the interior evaluation of vol um e integral 
computations. 

3.4 Artificial Dissipation 

Though the scheme is linearly stable, non-linear numerical instabilities could arise 
in strongly convective flows. Various artificial dissipation terms have been developed in 
the literature to suppress the numerical instabilities. The features we seek for artificial 
dissipation terms are that they only suppress numerical instabilities, they be small er than 
the real viscous terms, they are of higher order than the truncation terms and finally they 
should be implementable in the code without excessive computation. For this purpose, we 
choose the adaptive artificial dissipation terms of third order similar to those developed 
by Jameson^ 1 ) and others. These terms are included in the finite element code. A listing 
of the code is given in Appendix III. 


- 13 - 



4. Test Calculations 


4.1 Couette Flow 

The first test problem is the simulation of a one dimensional shear flow under pressure 
gradient. It has been computed with a uniform mesh of 2 x 6 x 2 linear (eight-node) 
elements with the following boundary conditions. 


u = t» = iu = 0at 
u = Uo, v = w = 0 at 
w = 0 at z = 0 and 
v = 0 at x = 0 and 


y = 0 plane 
y = 6 plane 
z = 2 plane 
x = 2 plane 


Q 

A favorable pressure gradient of — = — 1 is imposed. Fig. 4 shows the computed so- 

Ox 

lution with wail velocity Uq = 3. This problem has a simple exact solution as given in 
Schliching* 6 ). The computed solution agrees with the exact solution and the two are indis- 
tinguishable on the plot. For this simple problem, it takes very few time steps to reach a 
steady state solution starting from uniform flow conditions. The table of global and local 
correspondence of nodes, typical of finite element codes is also shown in Fig. 4. 

4.2 Laminar Boundary Layer Over a Flat Plate 

As another check case, laminar boundary layer over a flat plate has been computed 
with a stretched mesh of 4 x 6 x 1 linear elements. In this problem the convective terms are 
of the same order as some of the viscous terms. The finite element solution for a Reynolds 
number of Re = 10 4 , along with the boundary conditions and the mesh used are shown 
in Fig. 5. The computed solution agrees qualitatively with the exact solution even with a 
very coarse mesh. A converged solution can also be obtained for Re = 10 5 . 

4.3 Flow Over an Airfoil 

The locally implicit scheme for two dimensional Navier-Stokes equations with finite 
volume discretization is applied to compute airfoil flows. Calculations have been carried 
out with the code and comparisons have been made with experimental results. High 
Reynolds number viscous flows over an RAE 2822 airfoil have been computed from subsonic 
to transonic Mach numbers. An algebraic turbulence model developed by Baldwin and 
Lomax* 5 ) has been incorporated into the code. A body conforming C-grid (128 x 32) for 
an RAE 2822 airfoil is shown in Fig. 6. The mesh spacing normal to the airfoil is highly 
stretched to resolve turbulent viscous layer. The spacing ranges from .00005 to 3 chord 
lengths from inner to outer grid lines. Mach contours for turbulent flow at Mach number, 
M = 0.6, angle of attack, ot = 2.57 and Reynolds number, Re — 6.3 x 10 6 are shown in 
Fig. 7a. Fig. 7b shows the corresponding C p plot where numerical results are compared 
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with experimental values published by Cook, McDonald and Firmin^. The agreement of 
n um erical and experimental values for C p is reasonable for a relatively coarse gnd. Similar 
Mach contour and C p plots are presented for transonic flow case with M = 0.725, a = 2.92 
and Re = 6.5 x 10 6 in Figs. 8a and 8b. 

4.4 Flow in a Turn-around Duct 

As a test for the 3-D fini te element code, flow in an axisymmetric turnaround duct is 
computed at Mach number = 0.1. The schematic sketch of the turnaround duct is shown 
in Fig. 9. The geometry used corresponds to a test rig at Rockwell International which is 
shown in Fig. 10. A relatively coarse gnd of 8 x 15 x 2 elements are chosen. Since the 
flow is axisymmetric, 3 sectional planes with 2 elements in the circumferential direction 
are chosen and flow is set to be the same in each of the planes in the boundary conditions. 
The grid in one of the constant-angle planes and the computed velocity vectors are shown 
in Fig. 11 and a more detailed view of the velocity vectors in the bend region are shown 
in Fig. 12. The flow features are qualitatively correct. But a finer grid computation is 
necessary for quantitative comparisons with experimental results and it will be carried out 
later. 


- 15 - 


5. References 


1. Jameson, A., Baker, T. J., Weatherill, N. P., “Calculation of Inviscid Transonic Flow 
Over a Complete Aircraft”, AIAA-86-0103, AIAA 24th Aerospace Sciences Meeting, 
January 1986. 

2. Reddy, K. C., Nayani, S. N., “A Locally Implicit Scheme for Elliptic Partial Differential 
Equations”, presented at the SSME/CFD Working Group Meeting, NASA Marshall 
Space Flight Center, April 8-11, 1986. 

3. Reddy, K. C., Jacocks, J. L., “A Locally Implicit Scheme for the Euler Equations”, 
Proceedings of the AIAA 8th Computational Fluid Dynamics Conference, Honolulu, 
June 1987. 

4. Reddy, K. C., Reddy, J. N., Nayani, S. N., “Finite Element Solver for 3-D Compress- 
ible Viscous Flows”, Interim Report of Contract No. NASA8-36555, September 1987, 
submitted to NASA/MSFC, Marshall Space Flight Center, AL by The University of 
Tennessee Space Institute, Tullahoma, TN. 

5. Baldwin, B. S., Lomax, H., “Thin Layer Approximation and Algebraic Model for 
Seperated Turbulent Flows”, AIAA Paper 78-257, January 1978. 

6. Schlichting, H., Boundary Layer Theory , Pergamon Press, 1955. 

7. Reddy, J. N., An Introduction to the Finite Element Method , McGraw-Hill Book 
Company, 1984. 

8. Cook, P. H., McDonald, M. A., Firmin, M. C. P., “Airfoil RAE 2822 - Pressure Dis- 
tributions, and Boundary Layer and Wake Measurements”, AGARD-AR-138, 1979. 


- 16 - 






0 R= C 

x C/16 

A C/32 

1 Cl 64 

* C/2000 








































Fig. 6 Computational Grid for Viscous Flows 
RAE 2822 Airfoil - C grid (128 x 32) 
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Fig. 7a Mach Number Contours for Viscous Flow 

RAE 2822 Airfoil -M 00 = 0.6, a = 2.57°, Re = 6.3 x 10 6 
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Fig. 7b Numerical and Experimental Pressure Coefficients 

RAE 2822 Airfoil - Af* = 0.6, o = 2.57°, Re = 6.3 x 10 6 


- 24 - 



1 


0 . 


- 1 . 



Fig. 8a Mach Number Contours for Viscous Flow 

RAE 2822 Airfoil - = 0.7.25, a = 2.92°, Re = 6.5 x 10 6 
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Fig. 8b Numerical and Experimental Pressure Coefficients 

RAE 2822 Airfoil -M ao = 0.725, a = 2.92°, Re = 6.5 x 10 6 
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HOT-W I RE-PROBE PORTS LOCATED TO ALLOW 
DETECTION OF FLOW CHANGES THROUGH BEND 
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Fig. 10 Geometry of a Test Rig for a Turnaround Duct 


- 28 - 


TRAVERSE LINES FOR OTHER LOCATIONS SHOWN. 
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APPENDIX I 


The details of the Unsteady Compressible Navier— Stokes equations, which are used in 
the finite element code are given below. The equations are written in conservation form as 


where 


{U} = 


r p ' 

pu 

pv 

II 

i } 

»g + 

f S. 

II 

pw 


{-L’V + g) 

{ i \pe + p) ) 

, pe j 





£ — VjT, Tij — ~~ 2/iejj 


P = (7 - 1) pe - | (u 2 + v 2 + TP 2 )] tij = i (Uij + Vj t i) 

The viscous and inviscid fluxes are given by 


F v 

(5 x 3) 


' 0 

0 

0 ' 


' pu 

pv 

pw \ 

Til 

Tl2 

Tl3 

> , F 1 = . 

pu 2 +p 

puv 

puw 1 

T21 

T 22 

T23 

pvu 

pv 2 +p 

pvw j 

T31 

T32 

T33 

(5x3) 

pwu 

pwv 

pw 2 + p I 

D\ 

D 2 

D3 j 

- u(/>£ + p) 

v(pe + p) 

w(pe + p) J 


P= (7 ~ 1) c- ^(u 2 +u 2 +u> 2 )] ( p = pRT ), e = pe 


• Sutherland’s theory of viscosity : 


p = po 



( Tp + S\ \ 

Vr + sJ 


Si = constant (= 110 °K for air) 

• Properties of air at 20 °C(= T 0 ) and atmospheric pressure (pi = 1 atm) 
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Po = 1.205 Kg/m 3 

po = 0.101325 x 10 6 iV/m 2 

T 0 = 20 V = 293 °K 

sfhc) 

PO = 17.9 x 10" 6 (Pa - Sec) 
k = 2.5 x lO~ 2 (W/m - °K) 

P r = 0.72 
a = 0.208 
7 = 1.402 


AUXILIARY RELATIONS 


p = Pressure (iV/m 2 ) 

T = Temperature ( °K) 



C p = Specific heat at constant pressure 
C v = Specific heat at constant volume 
R = Gas constant ( N • m/Kg — °K) 
k = Thermal conductivity (W/m — °K ) 
po = Reference viscosity (Pa — Sec.) 

T 0 = Reference temperature ( °K) 
po = Reference density (Kg/m 3 ) 
p = pRT 



C v = 


R 

7-1 


a = Thermal diffusitivity, = — — 

pC p 

P r = Prandtl number = 

k 

Mgo — Mach number = 
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APPENDIX II 

Details of Finite Element Equations 

The details of fini te element equations which approximate the Navier— Stokes equations 
are given below. In equation (3.8) the residual {TV } has two parts. One is a volume 
integral, T Z v and the other is a surface integral, 71,. 

{TV} = {R v } + {K,} 


where 


(k.) = / (*] 7 '{F n )as 

Jd(l' 


The components of {7^„} for which corresponds to a node I are given by 


where 
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For defining the components of {72.,} we write 


F n dS = F • ndS = F -dS 


ft ( d(y,z) d{z,x) d(x,y)\ 

~ F '{ w^rW^rWv)) didri 
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as derived in equation (11) of the last report^ , for a typical surface, say C 1 of an 
element. 


Denote 


(V„V 2 ,V 3 )= ( 


d(y,z) 9(z,x) d(x,y) 

W.vr mv) 


) 


Now the components of {7?., } for which corresponds to a node /, for a typical surface 

£ = 1 of an element can be written as 


K = 


<f {V X U 2 + V 2 U 3 + V 3 U 4 )*rdt dr) 
Jan • 






^jdi dr] 



- 35 - 



where 


« ■ jL {&'"• * m * ft m m * ft m + ■ ' M 
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a^ 7 dQ d*idQ a^jag 
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V V x 
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1) 

p?)] 


$lid£ dr} 


Components of {7^j } for other surfaces of an element can be written similarly. 

The coefficient C of equation (3.13) has volume integrals of the derivatives of viscous flux 
terms. The details of those integrals are given below. 


Denote 

t dP° ^ 

/ «*?«» • s m— dir = N °m,i 

JCl* UU <*,J 

Subscript {ND) corresponds to the local index t of the global node ND in element e. These 
integrals can be written as 

^i = 0 

ft «£i ± (!i\ . ?hjL (!i\ + riifl dv 

[3 dx dz\vj dy dy \u,J dz dz 

± f*A - JL [«*i _ »..«!■ JLl etc 
dx{u,J~u,[dx ‘ &x u t \' ’ 



where 
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ooooooonooooooooooooooooooonoooooooooooooooooooooooooooooooooooooo 


FILE NAME: COMPR3D VERS 3; DATE: FEB. 22, 1988; LINES: 2066 


I I 
I FINITE-ELEMENT ANALYSIS OF FLOWS OF VISCOUS, COMPRESSIBLE I 
I FLUIDS IN THREE-DIMENSIONAL ENCLOSURES. I 
I I 


THIS PROGRAM IS DEVELOPED BY PROFESSORS J. N. REDDY OF 
VIRGINIA POLYTECHNIC INSTITUTE AND K. C. REDDY OF THE 
UNIVERSITY OF TENNESSEE SPACE INSTITUTE. THE PROGRAM IS 
UNDER CONTINUOUS DEVELOPMENT DURING APRIL '86 TO PRESENT. 
UNAUTHORIZED USE OF THE PROGRAM IS PROHIBITED. 

DEVELOPED: APRIL 1986 - PRESENT 


DESCRIPTION OF THE VARIABLES 


CFL THE COURANT-FRIEDRICHS-LEVY NUMBER 

ELXYZ ARRAY OF ELEMENT COORDINATES OF NODES 

IBNDC ARRAY OF BOUNDARY NODES FOR DIFFERENT 

VARIABLES 


I ORDER ORDER OF THE EQUATIONS TO BE SOLVED 

ISTART RESTART INDEX (1-RESTART; 0-NEW START) 

KELSUR A TWO-DIMENSIONAL ARRAY THAT CONTAINS ELEMENT 

NUMBER AND LOCAL NUMBER OF ITS SURFACE THAT 
REQUIRES FLUX COMPUTATION: 

KELSUR ( I, 1) -GLOBAL ELEMENT NUMBER OF THE 
GLOBAL I-TH SURFACE 

KELSUR (I, 2) -LOCAL SURFACE NUMBER OF THE 
GLOBAL I-TH SURFACE 

KNDSUR A TWO-DIMENSIONAL (M BY 4) ARRAY WHICH CONTAINS 

GLOBAL SURFACE NUMBERS SURROUNDING A NODE THAT 
REQUIRES FLUX COMPUTATION. HERE M DENOTES THE 
NUMBER OF NODES REQUIRING FLUX COMPUTATION: 

KNDSUR (I, J) -GLOBAL NUMBER OF THE LOCAL J-TH 
SURFACE ASSOCIATED WITH THE I-TH 
BOUNDARY NODE THAT REQUIRES FLUX 
COMPUTATION. 


MEN MAXIMUM NUMBER OF ELEMENTS AT A NODE 

MNE MAXIMUM NUMBER OF NODES PER ELEMENT 

NDF NO. OF UNKNOWNS AT EACH NODE 


NDSURF ARRAY CONTAINING THE SEQUENTIAL NUMBER OF THE 

BOUNDARY NODES WHICH REQUIRE FLUX COMPUTATION 
OR CONTAINING ZERO: 

NDSURF (I) -0, IF NO SURFACES AROUND THE I-TH 
NODE REQUIRES FLUX COMPUTATION. 

NDSURF (I) -J, IF THE I-TH NODE REQUIRES FLUX 
COMPUTATION; HERE J DENOTES THE SEQUENTIAL 
NUMBER OF NODE I IN THE LIST OF SURFACES THAT 
REQUIRE FLUX COMPUTATION. 

NELEM CONNECTIVITY MATRIX RELATING GLOBAL NODE TO 



oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


ELEMENTS AROUND THE NODE: 

NELEM ( I, M) -GLOBAL ELEMENT NUMBER CORRESPONDING 
TO THE M-TH LOCAL ELEMENT SURROUNDING GLOBAL 
NODE I (MAXIMUM VALUE OF M IS 8) . 


NEM NUMBER OF ELEMENTS IN THE MESH 

NGP NUMBER OF GAUSSIAN POINTS 

NMSH INDICATOR FOR GENERATING MESH: 


NMSH-0, MESH INFORMATION IS TO BE READ 
NMSH>0, MESH IS GENERATED BY THE PROGRAM 

(ONLY FOR PRISMATIC AND TAD DOMAINS) 


NNM NUMBER OF NODES IN THE MESH 

NODES BOOLEAN MATRIX RELATING LOCAL NODES TO GLOBAL 


NODES OF ELEMENTS: 

NODES (N, J) -GLOBAL NODE NUMBER CORRESPONDING TO 
THE J-TH LOCAL NODE OF ELEMENT N. 


NSURF TOTAL NUMBER OF SURFACES THAT REQUIRE 

FLUX COMPUTATION 
NTMSTP NO. OF TIME STEPS 

U ARRAY OF FIVE PRIMARY UNKNOWNS: 

RHO, RHO*U, RHO*V, RHO*W, RHO*E 

X,Y, Z GLOBAL COORDINATES OF THE NODES 


SUBROUTINES USED 

BCUPDT UPDATES THE BOUNDARY CONDITIONS AT THE END OF 

EACH ITERATION OR TIME STEP. 

BNDRY GENERATES ARRAY ' KNDSUR' , CONTAINING SURFACES 

REQUIRING FLUX COMPUTATION. 

COEFNT GENERATES THE COEFFICIENT VALUES OF EACH 

VARIABLE AT EACH NODE OF THE MESH. 

DISPTN COMPUTES THE DISSIPATION MODEL. 

DSFSUR COMPUTES THE DERIVATIVES OF THE SHAPE FUNCTIONS 

AT GAUSS POINTS OF A SURFACE. 

FLUXES COMPUTES FLUX FOR EACH VARIABLE AT EACH NODE 

OF THE MESH. 

GCSURF GENERATES ARRAY 'GC' , WHICH CONTAINS THE 

DERIVATIVE OF X(I) W.R.T. XI (J) . 

GMETRY GENERATES ARRAYS 'SF' , 'CNST', 'GDSF' AND 'VOL' 

GLOBALLY. 

INTIAL GENERATES INITIAL CONDITIONS ON BOUNDARY FA CE S. 

INVDET COMPUTES THE INVERSE OF THE JACOBIAN MATRIX. 

MATMUL COMPUTES THE PRODUCT OF TWO MATRICES. 


oooooooo o o o o n o o r> o o o o o o o o o o o o o o 


SHAPEL EVALUATES THE SHAPE FUNCTIONS AND THEIR DERIVA- 

TIVES AT THE GUASS POINTS. 

SURFGM COMPUTES COMPONENTS OF THE UNIT NORMAL AT 

GAUSS POINTS OF EACH BOUNDARY SURFACE. 

TADMSH GENERATES THE MESH ( X, Y AND Z COORDINATES AND 

ARRAY 'NODES') FOR THE TURN-AROUND-DUCT (TAD). 


IMPLICIT REAL* 8 <A-H,0-Z) 

PARAMETER (NNM-432 , NEM-240, MXE-8 , NGP-2 , NDIM-3, NPE-8, NDF-5, 

1 NBS-600) 

DIMENSION X(NNM) , Y (NNM) , Z (NNM) , TITLE (20) , UOLD (NNM, € ) ,U(NNM, 6) , 

2 NODES (NEM, NPE) ,NE LEM (NNM, MXE) , ELXYZ (NPE, NDIM) , E0 (NNM) , 

3 I ORDER (NDF) , DIS4 (NNM, 6) , DC4 (NNM) ,DELU (NPE, 6) , AMU (NNM) , 

4 GDSF (MXE, NPE, NGP, NGP, NGP, NDIM) , GNORM (NDIM, NBS , NGP , NGP ) , 

5 SF (NPE, NGP, NGP, NGP) ,CNST (MXE, NGP, NGP, NGP) ,EMU (NPE) , 

6 VOLND (NNM) ,VOL (MXE) ,DSURF (NDIM, NPE, 6, NGP, NGP) , 

7 ELU (NPE, 6) , IEL (MXE) , IBNDC (NNM, NDF) ,MINDX(NPE) , 

8 KELSUR (NBS , 2 ) , KNDSUR (NBS , 4 ) , NDSURF (NNM) 

COMMON/GMT /SN2 2 (8, 8) , SN33 (8, 8) , SN44 (8, 8) , SN55 (8, 8) 
COMMON/DTA/GAMA, AMUO , TEMPO , SI, R0 , GPR, GAM1, CFL 
DATA IORDER/1, 2,3,4, 5/ 

DATA IN, IT/5, 6/ 


I I 

(PREPROCESSOR | 

I I 

READ (5,2000) TITLE 

READ ( 5 , * ) ISTART,NMSH, ITER, NTMSTP, CFL, RLXOUT,RLXIN 
READ (5, *) AMUO, TEMP 0,S1,R0, GAMA, PR, AMACH,DNST0 
IF (NMSH . EQ . 0 ) GOTO 5 


CALL TADMSH (X, Y, Z, IBNDC, KELSUR, NODES, NSURF, NNM, NBS, NDF, NEM, NPE) 


GOTO 10 

5 READ ( 5 , * ) ( (NODES (I, J) , J— 1, 8) , I— 1,NEM) 

READ (5, *) ( (NELEM(I, J) , J-1,MXE> ,1-1, NNM) 

READ (5, *) (X(I) ,Y(I) ,Z (I), 1-1, NNM) 

READ (5,*) ( (U(I, J) , J-1,NDF) ,1-1, NNM) 

READ (5,*) NSURF 
IF (NSURF.EQ.O)GOTO 10 

READ (5, *) ( (KELSUR(I, J) , J-1,2) ,1-1, NSURF) 

READ (5,*) ( (IBNDC (I, J) , J— 1, 5) , I-1,NNM) 

END OF THE INPUT DATA 


OPEN THE OUTPUT FILE IN WHICH THE DATA IS TO BE STORED. 

THE NAME OF THE FILE IS 'TEST' AND THE DATA IS STORED IN THE FORM 
OF BINARY NUMBERS. 

10 CONTINUE 
IREC-30000 

OPEN (UNIT-08, FILE-' TEST' , STATUS-' NEW' , ACCESS-' DIRECT' , 

# FORM-' UNFORMATTED' , RECL-IREC, ACTION-' READ WRITE' ) 

IF ( I START . EQ . 1 ) THEN 

OPEN (UNIT-07, FILE-' RSTART' , STATUS-' OLD' , ACCESS-' DIRECT' , 

# FORM-' UNFORMATTED' , RECL-IREC, ACTION-' READWRITE ' ) 



ooo o on o ooo ooo ooo 


END IF 

GENERATE ARRAY 'NELEM' USING ARRAY ' NODES' 

DO 40 1=1, NNM 
DO 15 L=1,MXE 
15 NELEM ( I , L) =0 
ICNT=0 

DO 30 J=l, NEM 
DO 20 K=l, 8 
JK=NODES ( J, K) 

IF ( I . NE . JK ) GOTO 20 
ICNT=ICNT+1 
NELEM(I, ICNT) =J 
IF ( ICNT . EQ . MXE ) GOTO 40 
GOTO 30 
20 CONTINUE 
30 CONTINUE 
40 CONTINUE 

DEFINE FIXED PARAMETERS 

NGP T=NGP *NGP *NGP 
GAM1=GAMA-1 . 0 
GPR=GAMA/PR 

INITIALIZE THE FLOW FIELD 
NINIT=0 

IF (I START .EQ. 0) THEN 


CALL INTIAL (NDF, NNM, AMACH, AMU0 , TEMPO , SI , R0 , GAMA, PR, U, DNST0 ) 


CALL BCUPDT (NNM, GAMA, R0, TEMPO, U,DNST0) 


ELSE 

READ (07,REC=1) NINIT,U 
END IF 

NTMSTP = NTMSTP + NINIT 
NIN I T=N INI T+ 1 
DO 50 11=1, 6 
DO 50 JJ=1 , NNM 
50 UOLD (JJ, II) =U ( JJ, II) 

WRITE OUT INPUT DATA 

WRITE (IT, 2600) 

WRITE (IT, 2500) 

WRITE (IT, 2600) 

WRITE (IT, 3000) TITLE 

WRITE ( IT, 2 1 0 0 ) AMU0 , TEMPO , SI , R0 , GAMA, PR, DNST0 
WRITE ( IT, 22 0 0 ) ITER, NTMSTP, CFL, RLXOUT, RLXIN 
WRITE (IT, 7 41) AMACH 

741 FORMAT (1 OX, 'FREE STREAM MACH NUMBER =',E10.4) 
WRITE(IT, 3500) 

DO 70 I = 1, NEM 

70 WRITE (IT, 4000) I, (NODES (I, J) , J=l, 8) 

WRITE (IT, 4500) 

DO 80 I = 1, NNM 

80 WRITE (IT, 4000) I, (NELEM(I, J) , J=1,MXE) 

WRITE (IT, 5500) 

DO 90 I = 1, NNM 

90 WRITE (IT, 5000) I,X (I) , Y (I) , Z (I) 

WRITEdT, 6100) 

DO 100 1=1, NNM 



100 WRITE (IT, 6500) I, (U(I, J) , J=l, 5) 

WRITE(IT, 6200) 

DO 110 1=1, NNM 

110 WRITE (IT, 4000) I, (IBNDC (I, J) , J=l, 5) 

WRITE(IT, 6300) 

WRITE (IT, 4000) ( (KELSUR (I, J) , J=l, 2) , 1=1, NSURF) 

C 

C FIND MAX. NO. OF NODES PER EACH ELEMENT, COMPUTE ELEMENTAL 

C VOLUMES, SHAPE FUNCTIONS AND THEIR GLOBAL DERIVATIVES, AND 

C THE PRODUCT OF THE WEIGHTS AND THE DETERMINANT OF THE JACOBIAN 

C MATRIX FOR EACH GAUSS POINT OF EACH ELEMENT. 

C 

DO 155 ND=1,NNM 
C 

C COMPUTE THE NUMBER OF ELEMENTS AROUND NODE ' ND' 

C 

DO 115 J=1,MXE 
IF (NELEM (ND, J) .EQ.O)GOTO 120 
115 CONTINUE 
J=MXE+1 
120 NUMEL=J-1 
C 

C INITIALIZE THE ARRAYS 

C 

VOLND (ND) =0 . 0 
DC 4 (ND) =7*NUMEL 
C 

C COMPUTE ARRAY ' IEL' WHICH CONTAINS LOCAL NODE CORR TO NODE ND 

C 

DO 150 N=1,NUMEL 
NEL=NELEM (ND, N) 

DO 140 1=1, NPE 
NI=NODES (NEL, I) 

IF(NI.EQ.ND) IEL (N) =1 
ELXYZ (1,1) =X (NI ) 

ELXYZ (1,2) =Y (NI ) 

140 ELXYZ (I, 3) =Z (NI) 


CALL GMETRY (NNM, NEM, MXE, N, NPE, NGP, ELXYZ, SF, GDSF, CNST, VOL, 

1 NDIM, IEL (N) ) 

C 

150 VOLND (ND) =VOLND (ND) +VOL (N) 

WRITE (08, REC=ND) ND, CNST, GDSF, VOL,NUMEL, IEL, SN22, SN33, 

1 SN44, SN55 

* PRINT*, ND, CNST (1,1, 1,1), GDSF (1, 1, 1, 1, 1, 1) , VOL(l) 

155 CONTINUE 

C* WRITE (IT, 8000) (VOL (I) ,1=1, NEM) 

C 

CALL BNDRY (NBS, NEM, NNM, NPE, NSURF, NODES, KELSUR, NDSURF, KNDSUR) 
CALL DSFSUR(DSURF, NGP, NPE, NDIM) 


* WRITE (IT, 1000) 

* WRITE (IT, 4000) ( (KELSUR ( I, J) ,J=1, 2), 1=1, NSURF) 

* WRITE (IT, 4000) (NDSURF (I) , 1=1, 16) 

* WRITE (IT, 4000) ( (KNDSUR (I, J) , J=l, 4) , 1=1, NSURF) 

C 

DO 180 NDS=1, NSURF 
KE=KELSUR (NDS, 1) 

K1=KELSUR (NDS, 2 ) 

DO 160 1=1, NPE 
NI=NODES (KE, I) 

ELXYZ (1,1 )=X(NI) 

ELXYZ (1,2) =Y (NI) 

160 ELXYZ (1,3) =Z (NI ) 

C 

180 CALL SURFGM(K1, NDS, ELXYZ, DSURF, GNORM, NBS, NGP, NPE, NDIM) 



c 

c 
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BEGIN THE DO-LOOP ON THE 

NUMBER 

OF 

TIME 

STEPS 

TO 

COMPUTE THE SOLN 


ERROR=0 . 0 

DO 800 ITMSTP-NINIT , NTMSTP 
WRITE (IT, 6000) ITMSTP 
DO 190 1=1, NNM 
TEMP=U (1,6) /RO/U (1,1) 

190 AMU (I) =AMU0* ( (TEMP/TEMPO) **1.5) * ( (TEMP0+S1) / (TEMP+S1) ) 

CALL SUBROUTINE 'DISPTN' TO COMPUTE GLOBAL ARTIFICIAL DISSIPATION 


CALL DISPTN (NNM, NEM, MXE, X, Y, 2, U, DC 4 , NODES, NELEM, DIS4 , NPE, 
* E0 , VOLND) 


SYMMETRIC NONLINEAR GAUSS-SEIDEL ITERATION LOOP BEGINS HERE 

ITMAX=2 * ITER 
DO 700 ITR=1, ITMAX 
IF (MOD ( ITR, 2 ) . EQ . 1 ) THEN 
NBEGIN=1 
NEND=NNM 
NINC-1 
ELSE 

NBEGIN=NNM 

NEND=1 

NINC=-1 

ENDIF 

WRITE (IT, 4007) ITR, ITMAX 

BEGIN THE DO-LOOP ON THE NUMBER OF NODES TO COMPUTE THE SOLUTION 

DO 600 ND=NBEGIN, NEND, NINC 
WRITE (IT, 4006)NBEGIN, NEND, NINC, ND 

COMPUTE THE NUMBER OF ELEMENTS (NUMEL) SURROUNDING A NODE 

READ (08, REC=ND) ID, CNST, GDSF, VOL, NUMEL, IEL, SN22, SN33, 

1 SN44, SN55 

IF (ID .NE .ND) THEN 

PRINT *, 'ERROR IN THE READ OF FILES' 

STOP 

ENDIF 

NSTART=1 

NLAST=5 

INCR-1 

DO 500 LOOP=l, 1 

DO-LOOP ON THE NUMBER OF CONSERVATION EQUATIONS BEGINS HERE 

DO 400 NEQ=NSTART, NLAST, INCR 

WRITE (IT, 4004) NSTART, NLAST, INCR,NEQ, LOOP 

LEQ=IORDER (NEQ) 

IF (IBNDC (ND, LEQ) . EQ.OJGOTO 400 

DO-LOOP ON NUMBER OF ELEMENTS SURROUNDING NODE 'ND' BEGINS HERE 
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GCM=0.0 
GCKVIS=0 .0 
GCKINV=0 . 0 
TCOEF=0 . 0 
TRES=0 . 0 
TFLX=0 . 0 

DO 300 N=l, NUMEL 
WRITE (IT, 4003) NUMEL, N 
NEL=NELEM(ND,N) 

TRANSFER GLOBAL INFORMATION TO ELEMENT 'NEL' 

DO 260 1=1, NPE 
MINDX ( I ) =0 
NI=NODES (NEL, I) 

EMU (I) =AMU (NI) 

IF (NINC . EQ . 1 .AND. NI . GE .ND) MINDX ( I) =1 
IF (NINC.EQ.-l .AND. NI . LE .ND) MINDX ( I ) =1 
DO 260 11=1, 6 

DELU ( I , II) =U (NI, II) -UOLD (NI, II) 

260 ELU (I, II) =U (NI, II) 

CALL SUBROUTINE ' COEFNT' TO COMPUTE THE COEFFICIENTS FOR THE EQN 


CALL COEFNT (IEL (N) , LEQ, N, NPE, NEM, NGP, ELU, SF, GDSF, CNST, VOL, RES, 
* CM, EMU, DELU, MINDX, CKINV, NDF, NDIM, NGPT, MXE) 


GOTO (271,272,273,274,275) , LEQ 

271 GCKVIS=0 . 0 
GOTO 276 

272 DO 282 J1=1,NPE 

282 GCKVIS=GCKVIS+SN22 (N, Jl) *MINDX(J1) 
GOTO 276 

273 DO 283 J1=1,NPE 

283 GCKVIS=GCKVIS+SN33 (N, Jl) *MINDX(J1) 
GOTO 276 

274 DO 284 J1=1,NPE 

284 GCKVI S =GCKVI S +SN 4 4 (N, Jl) *MINDX(J1) 
GOTO 276 

275 DO 285 J1=1,NPE 

285 GCKVIS=GCKVIS+SN55 (N, Jl) *MINDX(J1) 

276 CONTINUE 
GCM=GCM+CM 
GCKINV=GCKINV+CKINV 

300 TRES=TRES+RES 

GCKINV=GCKINV*8 . 0 /NUMEL 
GCKVI S=GCKVIS* AMU (ND) /U(ND,1) 

IF ( LEQ . EQ . 5 ) GCKVI S=GCKVIS*GPR 
TCOEF=GCM+DABS (GCKINV) +GCKVIS 
TCOEF=TCOEF +DC4 (ND) 

IF (NDSURF (ND ) . EQ . 0 ) GOTO 340 

DO 335 J=l, 4 

KG1=KNDSUR (NDSURF (ND) , J) 

IF (KG1 . EQ . 0 ) GOTO 340 
K1=KELSUR(KG1,2) 

KL=KELSUR (KG1, 1) 

DO 310 11=1, NPE 

IF (NELEM (ND, II) . EQ.KL) THEN 

NI=II 

GOTO 315 

ENDIF 

310 CONTINUE 

315 DO 330 11=1, NPE 

EMU (II) =AMU (NODES (KL, II) ) 

DO 320 Jl=l , NDF 
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320 ELU ( I 1 , Jl ) =U (NODES (KL, I 1 ) , J1 ) 
330 IF (NODES (KL, II) ,EQ.ND)LI=I1 


CALL FLUXES (LI , LEQ, NI , NPE, NGP, ELU, SF, GDSF, GNORM, K1 , KG1 , FLX, 

1 EMU,MXE,NBS,NDF,NDIM) 

335 TFLX=TFLX+FLX 
340 CONTINUE 

IF (LEQ. NE. 2) GOTO 350 
ERRORO =ERROR 

ERROR=DMAX 1 (ERRORO, DABS (TRES+TFLX) ) 

IF (ERROR . GT . ERRORO ) MAXND=ND 
350 CONTINUE 
C DIS4 (ND, LEQ) =0 . 0 

DU=- (TRES+TFLX-DIS4 (ND, LEQ) ) /TCOEF 
U (ND, LEQ) =U (ND, LEQ) +DU*RLXIN 

U (ND, 6) =GAM1* (U(ND, 5) -0 .5* (U (ND,2) *U (ND, 2) +U (ND, 3) *U(ND, 3) + 

* U (ND, 4 ) *U (ND, 4 ) ) /U (ND, 1 ) ) 

WRITE ( IT, 7500) LEQ, ND, TRES, TFLX, TCOEF, U (ND, LEQ) 

400 CONTINUE 

NTEMP=N S TART 
NSTART=NLAST 
NLAST=NTEMP 
INCR=- 1 * INCR 
500 CONTINUE 

* WRITE(6, 9999) ND, (U (ND, LI) , LI=1, 6) 

*9999 FORMAT (15, 6E15 . 7) 

600 CONTINUE 

END OF THE COMPUTATION FOR ALL NODES IN THE SWEEP 

NTEMP=NBEGIN 
NBEG IN=NEND 
NEND=NTEMP 
NINC=-1*NINC 

RESET THE VALUES AT INFLOW, OUTFLOW AND RADIAL SYMMETRY PLANES 


CALL BCUPDT (NNM, GAMA, R0, TEMPO, U, DNST0) 


700 CONTINUE 

RELAXATION OF THE UPDATED SOLUTION AND COMPUTATION OF PRESSURE 

DO 720 11=1,5 
DO 720 JJ=1, NNM 

U( JJ, II) =UOLD ( JJ, II) +RLXOUT* (U ( JJ, II) -UOLD ( JJ, II) ) 

720 UOLD (JJ, II) =U (JJ, II) 

DO 730 Jl=l, NNM 

U( Jl, 6) =GAM1* (U ( Jl, 5) -0 . 5* (U (Jl, 2) *U (Jl, 2) +U (Jl, 3) *U (Jl, 3) + 

* U( Jl, 4) *U( Jl, 4) ) /U (Jl, 1) ) 

730 UOLD (Jl, 6) =U (Jl, 6) 

C* WRITE (IT, 7000) ERROR, MAXND 

DO 750 1=1, NNM 

750 WRITE(IT, 6500)1, (U(I, J) , J=l, 6) 

600 CONTINUE 

OPEN (UNIT=09, FILE=' ROLD' , STATUS='NEW' , ACCESS=' DIRECT' , 

* FORM=' UNFORMATTED' , RECL=IREC, ACTION=' READWRITE' ) 

WRITE ( 0 9, REC=1 ) NTMSTP, U 

C 

STOP 


FORMAT S 



1000 FORMAT 
2000 FORMAT 
2100 FORMAT 
2 

3 

4 

5 

6 

7 

8 

2200 FORMAT 
2 

3 

4 

5 

6 

2500 FORMAT 
2600 FORMAT 
3000 FORMAT 
3500 FORMAT 
★ 

4000 FORMAT 

4002 FORMAT 

4003 FORMAT 

4004 FORMAT 

4005 FORMAT 

4006 FORMAT 

4007 FORMAT 

4008 FORMAT 
4500 FORMAT 

* 

5000 FORMAT 
5500 FORMAT 
6000 FORMAT 
6100 FORMAT 
6200 FORMAT 
6300 FORMAT 


(5X, ' ARRAYS: KELSUR, NDSURF AND KNDSUR : ' , / ) 

(20A4) 

(/ , 2X, 'P ROBLEM DAT A:',/ 

/,5X, 'REFERENCE VISCOSITY (AMU0) =',E12.4 f 

/, 5X, 'REFERENCE TEMPERATURE (TEMPO) =',E12.4, 

/, 5X, ' SUTHERLANDS CONSTANT (SI) “',E12.4, 

/,5X, 'GAS CONSTANT (R0) =',E12.4, 

/,5X, 'RATIO OF SPECIFIC HEATS (GAMA) =',E12.4, 

/ , 5X, ' PRANDTL NUMBER (PR) «',E12.4, 

/,5X, 'REFERENCE DENSITY (DNST0) =',E12.4,/) 

(/, 2X,'P ARAMETERS OF APPROX. :',/, 

/,5X, 'NUMBER OF ITERATIONS PER TIME STEP =' , 13, 

/,5X, 'NUMBER OF TIME STEPS (NTMSTP) =',13, 

' "" ‘ =',E12.4, 

=' , E12 . 4, 

=' , E12 . 4, / ) 


/ , 5X, ' THE C F L NUMBER (CFL) 

/,5X, 'OUTER RELAXATION PARAMETER (RLXOUT) 
/,5X,' INNER RELAXATION PARAMETER (RLXIN) . 


UT FROM PROGRAM COMPR3D' , /) 


M A T R I X:',/, 


* COMPUTATION : ' , / ) 
6500 FORMAT (I5,6E12.4) 
7000 FORMAT 
7500 FORMAT 

* 

8000 FORMAT 
END 


(/, 15X, 'O U T P 
(80 ('-')) 

(1H1, 20A4) 

(/ , 2X, ' C ONNECTIVITY 
2X, ' (ELEMENT-TO-NODES) ' , /) 

(15, 2X, 1115) 

200 
300 
400 
500 
600 
700 
800 

(/ , 2X, ' C O N N E 

2X, ' (NODE-TO-ELEMENTS) ' , /) 

(15, 3 (2X, E12 . 4) ) 

(/ , 2X, ' (X, Y, Z ) -C OORDINATE 
(/ , 2X, ' T I M E STEP =',I5,/) 

(/ ,2X, 'INITIAL FIELD VALUES:',/) 

(/,2X, 'SPECIFIED NODAL QUANTITIES (=0, SPECIFIED):',/) 

(/, 2 X, 'ELEMENT NUMBERS AND THEIR SURFACES THAT REQUIRE FLUX 

m -r xt _ / / v 


(5X, 'DO- LOOP 
(5X, 'DO- LOOP 
(5X, 'DO-LOOP 
(5X, 'DO-LOOP 
(5X, 'DO-LOOP 
(5X, 'DO-LOOP 
(5X, 'DO-LOOP 


915) 
915) 
915) 
915) 
915) 
:',/,9l5) 
:',/,9 15) 
C T I V I 


T Y 


ARRAY 


,/, 


OF NODES:',/) 


(/,5X, 'MAX. ERROR »' , E12 . 4, / , 5X, ' NODE NUMBER -',15,/) 

( /, 5X, ' LEQ =' ,I2,2X, 'NODE =' , 14, 2X, 'RESIDUAL—' ,E12 . 4, 2X, 
' FLUX=' , E12 . 4 , 2X, ' TCOEF=' , E12 . 4, 2X, ' SOLN . =' , E12 . 4 ) 

( 5X, ' VOLUME OF EACH ELEMENT : ' , / , 5X, 6E12 . 4 ) 


SUBROUTINE BCUPDT (NNM, GAMA, R0, TEMPO , U, DNST0) 

C 

IMPLICIT REAL *8 (A-H,0-Z) 

COMMON /MSH / ARCANG, NX, NY, NZ , NX1 , NX2 , NX3 
DIMENSION U (NNM, 6) 

C 

C 

C DEFINE FIXED PARAMETERS 
C 

ANX=0 . 0 

ANY=DS IN ( 0 . 5 * ARCANG) 

ANZ=DCOS (0 . 5* ARCANG) 

GAM1—GAMA-1 . 0 
NXX-NX+1 
NYY=NY+1 
NZZ=NZ+1 
C 

C SET THE NORMAL VELOCITY TO ZERO AT THE MIDPLANE 

C 

DO 30 IX=1 , NXX 



DO 30 IY=1 , NYY 

ND= (IX-1) *NYY*NZZ+NYY+IY 

U (ND, 3) =U (ND, 3) * (1 . 0-ANY*ANY) -U (ND, 4) *ANY*ANZ 
U (ND, 4) =-U (ND, 3) *ANY*ANZ+U (ND, 4) * ( 1 . 0-ANZ*ANZ) 

U (ND, 5) =U (ND, 6) /GAM1+0 . 5* (U (ND, 2) *U (ND, 2) +U (ND, 3) *U(ND, 3) + 

* U (ND, 4) *U (ND, 4) ) /U (ND, 1) 

C 

C RESET THE VALUES ON PARALLEL PLANES TO THOSE ON THE MIDPLANE 

C 

ND1=ND-NYY 
ND2=ND+NYY 
U (ND1, 1 ) =U (ND, 1 ) 

U(ND1, 2) =U (ND, 2) 

U (ND1, 3 ) =U (ND, 3) *ANZ-U (ND, 4) *ANY 
U (ND1, 4) =U (ND, 3) *ANY+U (ND, 4) *ANZ 
U(ND1,5)=U(ND,5) 

U (ND1, 6) =U (ND, 6) 

U (ND2, 1 ) =U (ND, 1) 

U (ND2, 2 ) =U (ND, 2 ) 

U (ND2, 3) =U (ND, 3) *ANZ+U(ND, 4) *ANY 
U (ND2, 4 ) =-U (ND, 3) *ANY+U (ND, 4) *ANZ 
U (ND2, 5) =U (ND, 5) 

U (ND2, 6) =U (ND, 6) 

30 CONTINUE 
C 

C RESET THE VALUES AT OUTFLOW BOUNDARY 

C 

DO 40 IZ=1, NZZ 
DO 40 IY=1 , NYY 

ND = IY + ( IZ— 1 ) *NYY + NX*NYY*NZZ 
U (ND, 6) =DNST0*R0*TEMP0*0 . 98 

U (ND, 5) =U (ND, 6) /GAM1+0 . 5* (U (ND, 2) *U (ND, 2) +U (ND, 3) *U(ND, 3)+ 

* U(ND,4)*U(ND,4))/U(ND, 1) 

40 CONTINUE 
C 

C SET CONSTANT TEMPERATURE ON THE WALLS 

C 

DO 60 KD = 1, NX-1 
ND1 = (NYY*NZZ) *KD + 1 
DO 50 JZ ■= 1, NZZ 
ND = ND1 + (JZ-1) *NYY 
U (ND, 6) =U (ND, 5 ) *GAM1 
U (ND, 1 ) =U (ND, 6) / (R0*TEMP0) 

C 

NN - ND + NY 
U (NN, 6) =U (NN, 5) *GAM1 
U (NN, 1 ) =U (NN, 6 ) / (R0 *TEMP 0 ) 

50 CONTINUE 
C 

60 CONTINUE 
RETURN 
END 


SUBROUTINE BNDRY (NBS , NEM, NNM, NPE, NSURF, NODES, KELSUR, NDSURF, KNDSUR) 
IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION NODES (NEM, NPE) , KELSUR (NBS, 2) , KNDSUR (NBS, 4) , NDSURF (NNM) , 

* K (4) 

NCOUNT=0 
DO 10 1=1, NNM 
10 NDSURF (I)=0 
DO 20 L=l, 4 
DO 20 J=l, NSURF 
20 KNDSUR (J,L)=0 


C 
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DO 150 1=1 , NSURF 
KEL=KELSUR (1,1) 
KSRF=KELSUR(I,2) 

GOTO (30, 40, 50, 60,70, 80) , KSRF 

30 K(l)=NODES (KEL, 1) 

K (2) =NODES (KEL, 4) 

K ( 3 ) “NODES (KEL, 8) 

K(4) “NODES (KEL, 5) 

GOTO 90 

40 K(l) “NODES (KEL, 2) 

K (2) “NODES (KEL, 3) 

K (3) “NODES (KEL, 7) 

K (4) “NODES (KEL, 6) 

GOTO 90 

50 K (1) “NODES (KEL, 1) 

K (2) “NODES (KEL, 5) 

K (3) “NODES (KEL, 6) 

K (4) “NODES (KEL, 2) 

GOTO 90 

60 K(l) “NODES (KEL, 4) 

K (2) “NODES (KEL, 8) 

K(3) “NODES (KEL, 7) 

K(4) “NODES (KEL, 3) 

GOTO 90 

70 K(l) “NODES (KEL, 1) 

K (2) “NODES (KEL, 2) 

K (3) “NODES (KEL, 3) 

K (4) “NODES (KEL, 4) 

GOTO 90 

80 K(l) “NODES (KEL, 5) 

K (2) “NODES (KEL, 6) 

K (3) “NODES (KEL, 7) 

K (4) “NODES (KEL, 8) 

90 CONTINUE 

DO 120 J=l, 4 

IF (NDSURF (K ( J) ) . EQ . 0 ) THEN 
NCOUNT=NCOUNT+ 1 
NDSURF (K ( J) ) “NCOUNT 
KNDSUR (NCOUNT, 1 ) “I 
ELSE 

NC=NDSURF (K ( J) ) 

DO 100 JJ=2, 4 

IF (KNDSUR (NC,JJ) .EQ . 0) THEN 
KNDSUR (NC,JJ) =1 
GOTO 110 
ENDIF 

100 CONTINUE 
110 CONTINUE 
ENDIF 

120 CONTINUE 
150 CONTINUE 
RETURN 
END 


SUBROUTINE COEFNT ( IEL, LEQ, N, NPE, NEM, NGP , ELU, SF, GDSF, CNST, VOL, RES , 

CM, EMU, DELU, MINDX, CKINV, NDF, NDIM, NGPT, MXE) 


ELU (I, J) ELEMENT SOLUTION VECTOR (J-TH COMPO. AT I-TH NODE) 

SF (I, . . . ) SHAPE FUNCTION ASSOCIATED WITH THE I-TH NODE 

GDSF (N, J, .. I) .GLOBAL DERIVATIVE OF J-TH SHAPE FUNCTION 
WITH RESPECT TO X(I) COORDINATE 
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THIS IS A VECTORIZED VERSION OF THE SUBROUTINE COEFNT 
IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION SF (NPE, NGP, NGP, NGP) , CNST (MXE, NGP, NGP, NGP) f VOL(MXE) , 

2 GDSF (MXE, NPE, NGP, NGP, NGP, NDIM) , ELU (NPE, 6) , EMU (NPE) , 

3 U(6, 8) ,DU(7,3, 8) ,DU1 (7,3, 8) ,U1(6,8) , DELU (NPE, 6) f 

4 III (8) , JJJ(8) ,KKK(8) ,F(8, 8) ,DF<9,9,3) ,MINDX(NPE) , 

5 DQ1 (3) ,C (8) , GMU(8) 

COMMON/DTA/GAMA, AMUO , TEMPO , SI , RO , GPR, GAM1, CFL 

DATA 111/1,2,1,2,1,2,1,2/ 

DATA JJJ/1, 1,2, 2,1, 1,2,2/ 

DATA KKK/1, 1,1, 1,2, 2, 2, 2/ 


CM=0 . 0 
CK=0 . 0 
CKINV=0 . 0 
DLNGTH=0 . 0 
RES=0 . 0 
FMAS=0 . 0 

SPEED=DSQRT (ELU (IEL, 6) *GAMA/ELU (IEL, 1) ) 

DO 10 L=l, NGPT 

C (L) = CNST (N, III (L) , JJJ (L) , KKK (L) ) 

DO 10 1=1, NPE 

F (L, I) = SF (I, III (L) ,JJJ (L) , KKK (L) ) 

DF (L, I, 1) = GDSF (N,I,III(L), JJJ (L) , KKK (L) , 1 ) 

DF (L, 1,2) = GDSF (N, I, III (L) ,JJJ (L) , KKK (L) ,2) 

10 DF (L, I, 3) = GDSF (N, I, III (L) , JJJ (L) , KKK (L) , 3) 

TSPEED=SPEED+ (DABS (ELU (IEL, 2) ) +DABS (ELU (IEL, 3) ) +DABS (ELU (IEL, 4) ) ) / 
* ELU ( IEL, 1) 

DT=CFL* (VOL (N) ** (1 . /3 . ) ) /TSPEED 

EVALUATE THE SOLUTION AND ITS DERIVATIVES AT THE GAUSS POINT 

DO 40 J=l, NDF 
DO 40 L=1,NGPT 
SUM1=0 . 0 
SUM2=0 . 0 
SUM3=0 . 0 
SUM4=0 . 0 
DO 30 1=1, NPE 

SUM1=SUM1+DF (L, 1,1) *ELU(I, J) 

SUM2=SUM2+DF (L, 1,2) *ELU (I, J) 

SUM3=SUM3+DF (L, 1,3) *ELU (I, J) 

30 SUM4=SUM4+F (L, I) *ELU (I, J) 

DU ( J, 1, L) =SUM1 
DU ( J, 2, L) =SUM2 
DU ( J, 3 , L) =SUM3 
40 U ( J, L) =SUM4 
DO 50 J=2, 4 
DO 50 L=l, NGPT 
U1 ( J, L) =U ( J, L) /U (1, L) 

DU1 ( J, 1, L) = (DU (J, 1, L) -U1 (J, L) *DU (1, 1, L) ) 

DU1 ( J, 2, L) = (DU ( J, 2, L) -U1 (J,L) *DU(1,2,L) ) 

50 DU1 ( J, 3, L) = (DU ( J, 3, L) -U1 ( J, L) *DU (1, 3, L) ) 

COMPUTE MASS MATRIX TIMES DELU TERM 

DO 70 Jl=l, NPE 
DO 60 L=l, NGPT 
PROD=F (L, IEL) *F (L, Jl) *C (L) 

CM=CM+PROD*MINDX ( Jl ) 

6 0 FMAS =FMAS +P ROD *DELU ( Jl, LEQ ) 

70 CONTINUE 



c 

c 


COMPUTE INVISCID COEFFICIENT FOR INNER ITERATION- 

DO 90 L=1,NGPT 
CKINV=CKINV+ (DABS (DF (L, IEL, 1) * (DABS (U1 (2,L) ) +SPEED) ) 


1 
2 
3 

90 CONTINUE 


+ DABS (DF (L, IEL, 2) * (DABS (U1 (3,L) ) +SPEED) ) 

+ DABS (DF (L, IEL, 3) * (DABS (U1 (4,L) ) +SPEED) ) ) *C(L) 
*F (L, IEL) 


C 

C 

C 


COMPUTE RESIDUES ETC FOR A CONSERVATION EQUATION 
GOTO (100, 200, 300, 400, 500) , LEQ 


100 DO 110 L=1 , NGPT 

RES=RES- (DF (L, IEL, 1) *U (2, L) +DF (L, IEL, 2) *U(3, L) 

1 +DF (L, IEL, 3) *U (4, L) ) *C (L) 

110 CONTINUE 
GOTO 600 
C 

200 DO 240 L=l, NGPT 
SUM=0 . 0 

DO 220 1=1, NPE 
220 SUM=SUM+EMU (I) *F (L, I) 

240 <34U(L)=SUM 

DO 260 L=1 , NGPT 
U22=U (2, L) *U (2, L) 

U23=U (2, L) *U (3, L) 

U24=U (2, L) *U (4, L) 

U33=U (3, L) *U (3, L) 

U44=U (4, L) *U (4, L) 

PRES=GAM1* (U (5, L) -0 . 5* (U22+U33+U44) /U(1,L) ) 

AMU23=2 . 0*GMU (L) /3 . 0 
AMU43=2 . 0*AMU23 

RES=RES-C (L) * ( (U22+PRES*U ( 1 , L) +AMU23* (-2 . 0*DU1 (2, 1 . L) 

1 +DU1 (3, 2, L) +DU1 (4, 3, L) ) ) *DF (L, IEL, 1) 

2 + (U23-GMU (L) * (DU1 (3, 1, L) +DU1 (2,2,L) ) ) *DF(L,IEL,2) 

260 3 CONTINUE (U24_GMU (L) * (DU1 (4 ' lf L) +DU1 (2 ' 3 ' L > > > *DF (L, IEL, 3) ) /U (1, L) 
GOTO 600 

300 DO 340 L=l, NGPT 
SUM=0 . 0 

DO 320 1=1, NPE 
320 SUM=SUM+EMU (I) *F (L, I) 

340 GMU(L)=SUM 

DO 360 L=l, NGPT 
U22=U (2, L) *U (2, L) 

U23=U (2, L) *U (3, L) 

U33=U (3, L) *U (3, L) 

U34=U (3, L) *U (4, L) 

U44=U (4, L) *U (4, L) 

PRES=GAM1* (U (5, L) -0.5* (U22+U33+U44) /U(1,L) ) 

AMU23=2 . 0*GMU (L) /3 . 0 
AMU43=2.0*AMU23 

RES=RES-C (L) * ( (U33+PRES*U (1, L) +AMU23* (-2 . 0*DU1 (3,2 L) 

1 +DU1 (4, 3, L) +DU1 (2, 1, L) ) ) *DF (L, IEL, 2) 

2 + (U34-GMU (L) * (DU1 (4, 2, L) +DU1 (3, 3, L) ) ) *DF (L, IEL, 3) 

360 3 CONTINUE <U23 ***’ <L> * <DU1 (2, 2 ' L) +DU1 (3 ' 1 ' L) > > * DF < L » IE L, D ) /U (1, L) 
GOTO 600 


400 


420 

440 


DO 440 L=1 , NGPT 
SUM=0 . 0 

DO 420 1=1, NPE 
SUM=SUM+EMU (I) *F(L, I) 
GMU(L)=SUM 



DO 460 L=l, NGPT 
U22=U (2, L) *U (2, L) 

U24=U (2, L) *U (4, L) 

U33=U (3, L) *U (3, L) 

U34=U (3, L) *U (4, L) 

U44=U (4, L) *U (4, L) 

PRES=GAM1* (U (5, L) -0 . 5* (U22+U33+U44) /U(1,L) ) 

AMU2 3=2 . 0 *GMU ( L ) /3.0 
AMU43=2.0*AMU23 

RES=RES-C (L) * ( (U44+PRES*U (1, L) +AMU23* (-2 . 0*DU1 (4, 3, L) 

1 +DU1 (2, 1, L) +DU1 (3,2, L) ) ) *DF (L, IEL, 3) 

2 + (U24-GMU (L) * (DU1 (2, 3, L) +DU1 (4, 1, L) ) ) *DF (L, IEL, 1) 

3 + (U34-GMU (L) * (DU1 (3, 3,L)+DU1 (4,2,L) ) ) *DF (L, IEL, 2) ) /U(1,L) 
460 CONTINUE 

GOTO 600 


500 DO 540 L=1 , NGPT 
SUM=0 . 0 

DO 520 1=1, NPE 
520 SUM=SUM+EMU (I) *F (L, I) 

540 GMU(L)=SUM 

DO 560 L=1,NGPT 
U22=U (2, L) *U (2, L) 

U33=U (3, L) *U (3, L) 

U44=U (4, L) *U (4, L) 

PRES=GAM1* (U (5, L) -0 . 5* (U22+U33+U44) /U(1,L) ) 

AKH=GMU(L) *GPR 
AMU23=2 . 0*GMU (L) /3 . 0 
AMU43=2.0*AMU23 

DQ1 ( 1 ) =DU (5, 1, L) -U1 (2, L) *DU(2, 1,L) -U1 (3,L) *DU(3, 1,L) 

2 -U1 (4, L) *DU (4, 1, L) +DU (1, 1, L) * (-U (5, L) /U(l, L) 

3 +U1 (2, L) *U1 (2, L) +U1 (3, L) *U1 (3, L) +U1 (4, L) *U1 (4, L) ) 
DQ1 (2) =DU(5, 2, L) -U1 (2,L) *DU (2, 2, L) -U1 (3,L) *DU(3,2,L) 

2 -U1(4,L) *DU (4, 2, L) +DU (1, 2, L) * (-U (5, L) /U(l, L) 

3 +U1 (2, L) *U1 (2, L) +U1 (3, L) *U1 (3, L) +U1 (4, L) *U1 (4, L) ) 
DQ1 (3) =DU ( 5 , 3, L) -U1 (2, L) *DU (2, 3, L) -U1 (3, L) *DU(3, 3, L) 

2 -U1(4,L) *DU (4, 3, L) +DU (1, 3, L) * (-U (5, L) /U(1,L) 

3 +U1 (2, L) *U1 (2, L) +U1 (3, L) *U1 (3, L) +U1 (4, L) *U1 (4, L) ) 


REST 


2 

3 

4 

I 

2 

3 

4 


RES2 


RES 3 


2 

3 

4 


RES = RES 
560 CONTINUE 
600 CONTINUE 


= (U (2, L) * (U (5, L) +PRES) -AMU23*U1 (2, L) * (2 . 0*DU1 (2, 1, L) 
-DU1 (3, 2, L) -DU1 (4, 3, L) ) -GMU(L) * (U1 (3, L) * (DU1 (2,2, L) 
+DU1 (3, 1, L) ) +U1 (4, L) * (DU1 (2, 3, L) +DU1 (4, 1,L) ) ) 
-AKH*DQ1 (1) ) *DF (L, IEL, 1) 

= (U (3, L) * (U (5, L) +PRES) -AMU23*U1 (3, L) * (2 . 0*DU1 (3,2,L) 
-DU1 (4, 3, L) -DU1 (2, 1, L) ) -GMU(L) * (U1 (4, L) * (DU1 (3, 3, L) 
+DU1 (4, 2, L) ) +U1 (2, L) * (DU1 (3, 1, L) +DU1 (2,2, L) ) ) 
-AKH*DQ1 (2) ) *DF (L, IEL, 2) 

= (U (4, L) * (U(5, L) +PRES) -AMU23*U1 (4, L) * (2 . 0*DU1 (4, 3, L) 
-DU1 (2, 1, L) -DU1 (3,2, L) ) -GMU(L) * (U1 (2, L) * (DU1 (4, 1, L) 
+DU1 (2,3, L) ) +U1 (3, L) * (DU1 (4,2, L) +DU1 (3, 3,L) ) ) 
-AKH*DQ1 (3) ) *DF (L, IEL, 3) 


- (RES1+RES2+RES3) *C (L) /U (1, L) 


RES=RES+FMAS /DT 

CM=CM/DT 

RETURN 

END 


SUBROUTINE DISPTN (NNM, NEM, MXE, X, Y, 
* NPE, E0 , VOLND) 

IMPLICIT REAL* 8 (A-H,0-Z) 


Z,U,DC4, NODES, NELEM,DIS4, 



DIMENSION X (NNM) , Y (NNM) , Z (NNM) , U (NNM, 6) , NODES (NEM, 8) , EO (NNM) , 
2 NELEM<NNM,MXE) ,DIS4 (NNM, 6) , VOLND (NNM) ,DC4 (NNM) 

DATA KAPA2, KAPA4/0 .1,0.01/ 

DO 50 IE-1, 6 

DO 40 ND=1, NNM 

SUME0=0 . 0 

DO 20 NE=1,MXE 

IF (NELEM (ND, NE) . EQ.O)GOTO 30 

NEL=NELEM ( ND , NE ) 

DO 20 NP=1,NPE 
NI=NODES (NEL, NP) 

20 SUME0=SUME0+U (NI, IE) -U (ND, IE) 

NE-MXE+1 
30 CONTINUE 

DC 4 (ND) =7* (NE-1 ) 


40 DIS4 (ND, IE) =SUME0 
50 CONTINUE 

DO 60 ND=1, NNM 

DIS4 (ND, 5) =DIS4 (ND, 5) +DIS4 (ND, 6) 

60 DIS4 (ND, 6) =ABS (DIS4 (ND, 6) ) /U(ND, 6) *KAPA2 

COMPUTE THE FOURTH-ORDER DISSIPATION 


DO 150 IE-1,5 
DO 140 ND=1, NNM 
SUMDC-0 . 0 
E0 (ND) =0 . 0 
SUMD0=0 . 0 
ISW=1 


IF (DIS4 (ND, 6) .GT.KAPA4) ISW-0 
DO 120 NE=1,MXE 
NEL—NELEM (ND, NE) 

IF (NEL. EQ. 0) GOTO 130 
DO 100 NP=1 , NPE 
NI=NODES (NEL,NP) 

IF (NI . EQ . ND ) GOTO 100 
XL=X (NI ) -X (ND) 

YL=Y(NI)-Y(ND) 

ZL=Z (NI) -Z (ND) 


100 

120 

130 


EDGE -DSQRT (XL*XL+YL*YL+ZL*ZL) 

EPSLN= (VOLND (ND) +VOLND (NI) ) *0. 5/EDGE 

IF ( IE . EQ . 5) SUMDC=SUMDC+EPSLN* ( (DC4 (ND) —1 . 0) *KAPA4 

SUMD 0=SUMD 0 - (DIS4 (NI, IE) -DIS4 (ND, IE) ) *EPSLN*KAPA4 

CONTINUE 

CONTINUE 

CONTINUE 


*ISW+DIS4 (ND, 6) ) 
*ISW 


IF (IE .EQ .5) DC4 (ND) =SUMDC 
140 E0 (ND) =SUMD0 

DO 150 ND = 1, NNM 

150 DIS4 (ND, IE) =E0 (ND) +DIS4 (ND, IE) *DIS4 (ND, 6) 
RETURN 
END 


C 

c 

c 

c 

c 


SUBROUTINE DSFSUR (DSURF, NGP, NPE, NDIM) 


SU ^ TINE EVALUATES THE DERIVATIVES OF THE SHAPE FUNCTIONS 
AT THE GAUSS POINTS OF THE SURFACES OF AN ELEMENT 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION XNODE ( 8 , 3 ) , XYZ ( 3 ) , DSURF (NDIM, NPE 
DATA XNODE/ -1 . 0D0, 2*1 . 0D0, 2*-l . 0D0, 2*1 . 0D0 


6, NGP, NGP) , GAUSS (2) 

-1 . 0D0, 2*-l . 0D0, 2*1 . 0D0 



OOOOOOQOOO 


c 


1,2*-1.0DO,2*1.0DO,4*-1.0DO,4*1.0DO/ 


FCK (A, B, C) =0 . 125*A*B*C 
SQRT3=DSQRT (3 . 0D0) 

GAUSS (1)— 1 . 0D0/SQRT3 
GAUSS (2 ) =-GAUSS (1) 

DO 80 Kl=l, 6 
DO 60 NGPI=1,NGP 
DO 60 NGPK=1 , NGP 

GOTO (10, 10, 20, 20, 30, 30), K1 
10 XYZ < 1) = (-1 ) **K1 
XYZ (2) =GAUSS (NGP I) 

XYZ (3) =GAUSS (NGPK) 

GOTO 40 

20 XYZ (2) = (-1 ) **K1 
XYZ (3) =GAUSS (NGP I) 

XYZ (1)=GAUSS (NGPK) 

GOTO 40 

30 XYZ (3) = (-1) **K1 
XYZ (1)=GAUSS (NGP I) 

XYZ (2) =GAUSS (NGPK) 

40 DO 50 1=1, NPE 

XNP 1=XYZ (1) *XNODE (I, 1)+1 .0 
YNP 1=XYZ (2) * XNODE ( I , 2 ) +1 . 0 
ZNP 1=XYZ (3) *XNODE (I, 3) +1 .0 

DSURF ( 1 , I , K 1 , NGP I , NGPK ) =FCK ( XNODE (1,1), YNP 1 , ZNP 1 ) 
DSURF (2, I, Kl, NGP I, NGPK) =FCK (XNP 1, XNODE (1,2) , ZNP1) 
50 DSURF (3, I, K1,NGPI, NGPK) =FCK (XNP1, YNP 1, XNODE (1,3) ) 
60 CONTINUE 
80 CONTINUE 
RETURN 
END 


SUBROUTINE 

1 


FLUXES ( IEL, LEQ, N, NPE, NGP , ELU, SF, GDSF, GNORM, Kl, KG1 , FLX, 
EMU, MXE, NBS, NDF, NDIM) 


ELU (I, J) ELEMENT SOLUTION VECTOR (J-TH COMPO. AT I-TH NODE) 

SF(I, ...) SHAPE FUNCTION ASSOCIATED WITH THE I-TH NODE 

GDSF (N, J, .. I) .GLOBAL DERIVATIVE OF J-TH SHAPE FUNCTION 
^ T „ Trn , WITH RESPECT TO X(I) COORDINATE OF THE N-TH ELEMENT 

GDINT(I,J) ....INTERPOLATED GDSF ON SURFACE OF AN ELEMENT 
SFINT(I) INTERPOLATED SF ON SURFACE OF AN ELEMENT 


C 


C 

C 

C 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION SF (NPE, NGP, NGP, NGP ) , GDSF (MXE, NPE, NGP, NGP, NGP, NDIM) 

2 GD INT ( 8 , 3 ) , SF INT ( 8 ) , GNORM (ND IM, NBS , NGP , NGP ) , EMU ( NPE ) 

3 ELU (NPE, 6) ,DU (6, 3) , U ( 6) , U1 (6) , DU1 ( 6, 3) , DQ1 (3) , VECTR (3) 
COMMON/DTA/GAMA, AMU0, TEMPO , SI, R0 , GPR, GAM1, CFL 


K0= (Kl+1 ) / 2 
FLX=0 . 0 

SQRT3=DSQRT (3 . 0D0) 


DO-LOOP ON GAUSS INTEGRATION BEGINS HERE 

DO 200 JJ=1,NGP 
DO 200 KK= 1 , NGP 
AMU=0 . 0 


C 

C 

C 


EVALUATE THE COMPONENTS OF THE SURFACE NORMAL AT THE GAUSS POINTS 



IF (KO-2) 30, 40, 50 
30 NI=1 
NI1=2 
NJ=JJ 
NJ1=NJ 
NK=KK 
NK1=NK 
GOTO 60 
40 NJ=1 
NJ1=2 
NK=JJ 
NK1=NK 
NI=KK 
NI1=NI 
GOTO 60 
50 NK=1 
NK1=2 
NI=JJ 
NI1=NI 
NJ=KK 
NJ1=NJ 
C 

60 DO 70 1=1, NPE 

F1=SF(I,NI,NJ,NK) 

F2=SF (I, Nil, NJ1, NK1) 

SFINT ( I ) = ( (-1) **K1*SQRT3* (F2-F1) +F2+F1) /2.0 
F3=GDSF (N, I , NI , N J, NK, 1 ) 

F4=GDSF (N, I , Nil , NJ1 , NK1, 1) 

GDINT (I, 1) = ( (-1) **K1*SQRT3* (F4-F3) +F4+F3) /2 . 0 
F3=GDSF (N, I , NI , NJ, NK, 2 ) 

F4=GDSF (N, I , NI 1 , NJ1 , NK1, 2 ) 

GDINT (I, 2) = ( (-1) **K1*SQRT3* (F4-F3) +F4+F3) /2.0 
F3=GDSF (N, I , NI , NJ, NK, 3 ) 

F4=GDSF (N, I , NI 1 , NJ1 , NK1 , 3 ) 

GDINT {I, 3) = ( (-1) **K1*SQRT3* (F4-F3) +F4+F3) /2 . 0 
70 AMU=AMU+SFINT(I) *EMU(I) 

DO 100 J=l, NDF 
SUM1=0 . 0 
SUM2=0 . 0 
SUM3=0 . 0 
SUM4=0 . 0 
DO 80 1=1, NPE 

SUM1=SUM1+GDINT (1,1) *ELU (I, J) 

SUM2=SUM2+GDINT (1,2) *ELU ( I , J) 

SUM3=SUM3+GDINT (I, 3) *ELU (I, J) 

80 SUM4=SUM4+SFINT (I) *ELU (I, J) 

DU (J, 1) =SUM1 
DU (J, 2) =SUM2 
DU (J, 3) =SUM3 
100 U ( J) =SUM4 

U1 ( 2 ) =U ( 2 ) /U(l) 

U1(3)=U(3)/U(1) 

U1(4)=U(4)/U(1) 

DO 110 J=2, 4 

DU1 (J, 1) = (DU (J, 1)-U1 (J) *DU (1,1) ) 

DU1 ( J, 2) = (DU (J, 2) -U1 (J) *DU (1,2) ) 

110 DU1 ( J, 3) = (DU (J, 3) -U1 (J) *DU (1,3) ) 

VECTR(l) =GNORM (1, KG1, JJ, KK) 

VECTR (2) =GNORM (2, KG1, JJ, KK) 

VECTR ( 3 ) =GNORM ( 3 , KG1 , JJ, KK ) 

C 

C COMPUTE PRESSURE, TEMPERATURE, VISCOSITY USING THE SUTHERLAND'S 

C LAW, AND THE DIFFUSION CONSTANT AT THE GAUSS POINTS 

C 

U22=U (2) *U (2) 

U23=U (2 ) *U (3) 



c 

c 

c 


U24=U (2) *U(4) 

U33=U (3) *U (3) 

U34=U (3) *U (4) 

U44=U (4) *U (4) 

PRES=GAM1* (U (5) -0.5* (U22+U33+U44) /U(l) ) 

AKH=AMU*GPR 
AMU23=2 . 0*AMU/3 . 0 
AMU43=2 . 0*AMU23 

COMPUTE THE FLUX FOR EACH CONSERVATION EQUATION AT THE 


NODE 


140 


150 


GOTO (140, 150, 160, 170, 180), LEQ 

Swa**™ <U ( 2 ] * VECTR < 1 > +U ( 3 ) * VECTR ( 2 ) +U ( 4 ) * VECTR (3) ) *SFINT (IEL) 

FLX=FLX+ ( (U22+PRES*U (1) +AMU23* (-2 . 0*DU1 (2, 1) +DU1 (3, 2) +DU1 (4,3))) 

1 * VECTR ( 1 ) + (U23-AMU* (DU1 (3,1) +DU1 (2,2))) *VECTR(2) 

2 nrwn, onn ^ U24_AMU * < DU1 ’ D +DU1 (2,3))) *VECTR(3) ) *SFINT(IEL) /U (1) 
CiOTO Z U 0 

160 FLX=FLX+( (U33+PRES*U (1) +AMU23* (-2. 0*DU1 (3,2) +DU1 (4, 3)+DUl (2 1) ) ) 
1 *VECTR (2) + (U34-AMU* (DU1 (4,2) +DU1 (3, 3) ) ) *VECTR (3) 

2 GOTO 200 <U23_AMU * * DU1 2 ) +D U1 (3,1))) *VECTR(1) ) *SFINT(IEL) /U (1) 

170 FLX=FLX+ ( (U44+PRES *U (1 ) +AMU23* (-2. 0*DU1 (4, 3) +DU1 (2, 1)+DU1 (3 2) ) ) 

1 *VECTR(3) + (U24-AMU* (DU1 (2, 3) +DU1 (4, 1) ) ) *VECTR (1) 

GOTO 200 (U34_AMU * (DU1 (3 ' 3) +DU1 < 4 ' 2 ) ) ) *VECTR(2) ) *SFINT(IEL) /U(l) 

180 DQl (1) =DU(5, 1) -U1 (2) *DU(2, 1) -U1 (3) *DU(3, 1) -U1 (4) *DU(4, 1) 

^ +DU (1, 1) * (-U (5) /U (1) +U1 (2) *U1 (2) +U1 (3) *U1 (3) +U1 (4) *U1 (4) ) 

DQl (2) =DU (5, 2) -U1 (2) *DU (2,2) -U1 (3) *DU <3, 2) -ill (4) *DU <4, 2 

3 __ . V * * U ^ /U* 1 ) +U1 (2) *U1 (2) +U1 (3) *U1 (3) +U1 (4) *U1 (4) ) 

DQl (3) -DU ( 5 , 3) -U1 (2) *DU (2,3) -U1 (3) *DU(3, 3) -U1 (4) *DU(4, 3) 

2 +DU (1, 3) * (-U (5) /U (1) +U1 (2) *U1 (2) +U1 (3) *U1 (3) +U1 (4) *U1 (4) ) 
FLX-FLX+ ( (U (2) * (U(5) +PRES) -AMU23*U1 (2) * (2 . 0*DU1 (2, 1) -DU1 (3 2) 

2 -DUi (4,3) ) -AMU* (U1 (3) * (DU1 (2, 2) +DU1 (3, 1) ) +U1 (4) * (DU1 (2 3) 

3 +DU1 (4, 1) ) ) -AKH*DQ1 (1) ) *VECTR(1) ' 

4 + (U(3) * (U (5) +PRES ) -AMU2 3 *U1 (3) * (2 ,0*DU1 (3,2) -DUI (4, 3) 

5 -DUi (2, 1) )-AMU* (U1 (4) * (DUI (3,3) +DU1 (4,2) ) +U1 (2) * (DUI (3 1) 

6 +DU1 (2,2) ) ) -AKH*DQ1 (2) ) *VECTR(2) ' 

7 + (U (4) * (U (5) +PRES) -AMU23*U1 (4) * (2 . 0*DU1 (4,3) -DUI 12 11 

8 -DUI (3,2) ) -AMU* (U1 (2) * (DUI (4,1) +DU1 (2,3)) +U1 (3) * (DUI ( 4 2) 

200 CONTINUE^ 3 ' 3> 5 ) -A kh *DQK 3) ) *VECTR(3) ) *SFINT(IEL) /U(l) ' 

C* WRITE ( 6, 300) LEQ, FLX 

C 300 FORMAT (5X, 'LEQ =',I2,5X,'FLUX =',E12 4) 

RETURN 

END 


c SUBROUTINE GCSURF (GC, DSURF, ELXYZ, NGPI, NGPK, NGP, Kl, NDIM, NPE) 

S. GC(I,J) DERIVATIVE OF X(I) W.R.T. XI (J) 

C DSURF (I, J,K. .DERIVATIVE OF PSI ( J) W.R.T. XI(I), J=l, Np E 
C ON K-TH SURFACE OF MASTER ELEMENT 

IMPLICIT REAL*8 (A-H, O-Z) 

c DIMENSION ELXYZ (NPE, NDIM), DSURF (NDIM, NPE, 6, NGP, NGP), GC (NDIM, NDIM) 

DO 200 1=1, NDIM 
DO 200 K=l, NDIM 
SUM=0 . 0 

DO 100 J=1 , NPE 

100 SUM=SUM+DSURF (K, J,K1, NGPI, NGPK) *ELXYZ(J, I) 

200 GC ( I, K) =SUM 
RETURN 
END 



oooooooooooon 


SUBROUTINE GMETRY (NNM, NEM, MXE, N, NPE, NGP, ELXY2, SF, GDSF, CNST, VOL, 
1 NDIM, IEL) 


SF (I, II, JJ, KK) I-TH SHAPE FUNCTION AT THE (II, JJ, KK) -TH 

GAUSS POINT 

GDSF (N, I, II, JJ, KK, J) . .GLOBAL DERIVATIVE OF I-TH SHAPE FUNCTION 


WITH RESPECT TO THE X(J) COORDINATE 
FOR ELEMENT N 

DSF (I, J) LOCAL DERIVATIVE OF I-TH SHAPE FUNCTION 

_, TW „, T ,, WITH RESPECT TO J-TH LOCAL COORDINATE 

ELXYZ (I, J) J-TH GLOBAL COORDINATE OF I-TH NODE 

XYZdD I I-TH GAUSSIAN POINT 


C 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION SF (NPE, NGP, NGP, NGP) , CNST (MXE, NGP, NGP, NGP) ,VOL(MXE) 

2 GDSF (MXE, NPE, NGP, NGP, NGP, NDIM) , ELXYZ (NPE, NDIM) ,WT (2) 

3 GAUSS (2) , GJ (3, 3) , XYZ (3) , GJINV (3,3) ,DSF (3, 8) , GDSFL (3, 
SF L ( o ) 

COMMON / GMT /SN2 2 (8,8) ,SN33 (8, 8) , SN44 (8, 8) ,SN55 (8, 8) 

DATA NCOUNT/O/ 


8 ) 


C 

C 

c 


c 


c 


SQRT3=DSQRT (3 . 0D0 ) 
GAUSS (1) — 1 . 0D0/SQRT3 
GAUSS (2) =-GAUSS (1) 

WT (1) =1 . 0D0 
WT (2) =1 . 0D0 


DO-LOOP ON GAUSS INTEGRATION BEGINS HERE 


VOL (N) =0 . 0 


DO 50 J=l, NPE 
SN22 (N, J) =0 . 0 
SN33 (N, J) =0 . 0 
SN44 (N, J) =0 . 0 
50 SN55 (N, J) =0 . 0 
DO 200 11=1, NGP 
DO 200 JJ=1, NGP 
DO 200 KK=1, NGP 
XYZ (1) =GAUSS (II) 
XYZ (2) =GAUSS (JJ) 
XYZ (3) =GAUSS (KK) 






CALL SHAPEL (XYZ, SFL, DSF, NDIM, NPE) 

CALL MATMUL (DSF , ELXYZ, GJ, NDIM, NPE, NDIM) 

CALL INVDET(GJ, GJINV, DET) 

CALL MATMUL (GJINV, DSF, GDSFL, NDIM, NDIM, NPE) 

***************************** i 'i' i , i , i ,* i ' i '* i!i ' i ' i ' i ' i ' i!ili ' ic 

CNST (N, II, JJ,KK) =DET*WT (II) *WT (JJ) *WT (KK) 

DO 150 1=1, NPE 


******** 


★ * 


** 


****** 


****** 


SN22 (N, I)=SN22 (N, I) + (4 . 0/3 . 0*GDSFL (1, IEL) *GDSFL(1, I) + 

1 GDSFL (2, IEL) *GDSFL(2, I) +GDSFL (3, IEL) * GDSFL (3, I) ) * 

2 CNST (N, II, JJ, KK) 

SN33 (N, I) =SN33 (N, I) + (4 . 0/3 . 0*GDSFL (2, IEL) * GDSFL (2, I) + 

1 GDSFL (3, IEL) * GDSFL (3, 1) +GDSFL (1, IEL) * GDSFL (1, 1) ) * 

2 CNST (N, II, JJ, KK) ' ' 

SN44 (N, I) =SN44 (N, I) + (4 . 0/3 . 0 *GDSFL (3, IEL) * GDSFL (3 I) + 

1 GDSFL (1, IEL) *GDSFL (1, 1)+GDSFL (2, IEL) *GDSFL(2, I) ) * 

2 CNST (N, II, JJ, KK) ' " 

SN55 (N, I) =SN55 (N, I) + (GDSFL (1, IEL) *GDSFL(1, 1) + 

1 GDSFL (2, IEL) *GDSFL (2, I) +GDSFL (3, IEL) * GDSFL (3,1))* 

2 CNST (N, II, JJ, KK) ' " 



IF (NCOUNT . GT . 0 ) GOTO 100 
SF (I,II,JJ, KK) =SFL (I) 

100 GDSF (N, I, II, JJ,KK, 1) =GDSFL (1, I) 
GDSF (N, I, II, JJ, KK, 2) =GDSFL (2, I) 
150 GDSF (N, I, II, JJ,KK,3)=GDSFL(3,I) 
VOL (N) =VOL (N) +CNST (N, II, JJ, KK) 
200 CONTINUE 
NCOUNT =1 
RETURN 
END 


SUBROUTINE INTIAL (NDF, NNM, AMACH, AMU0 , TEMPO, SI, R0, GAMA, PR, U, DNST0 ) 

C INITIAL CONDITIONS FOR THE TURN-AROUND -DUCT PROBLEM 

C 

IMPLICIT REAL* 8 (A-H,0-Z) 

COMMON/MSH/ARCANG, NX, NY, NZ , NX1 , NX2 , NX3 
DIMENSION U (NNM, 6) 

C 

c 

C DEFINE FIXED PARAMETERS 

C 

GAM1=GAMA-1 . 0 
NYY=NY+1 
NZZ=NZ+1 
C 

C INITIALIZE THE FLOW FIELD 

C 

DO 10 J=2, 4 
DO 10 1=1, NNM 
10 U(I,J)=0.0 
C 

DO 20 IZ=1, NZZ 
DO 20 IY=2, NY 
ND = IY + (IZ-1) *NYY 
U (ND, 2) =-DSQRT (GAMA*R0*TEMP0) * AMACH 
IF (IY.EQ .2 .OR. IY .EQ. 8) U (ND, 2) =U (ND, 2) *0.1885 
IF (IY.EQ.3 .OR. IY.EQ.7) U (ND, 2)=U (ND,2) *0.5066 
IF (IY . EQ . 4 . OR. IY .EQ . 6) U (ND, 2 ) =U (ND, 2 ) *0.8393 
20 CONTINUE 
C 

C INITIALIZE THE MID PLANE 

C 

DO 30 IX = 2,NX1+1 
DO 30 IY = 2, NY 
ND = (IX-1) *NYY*NZZ+NYY+IY 
NDI= NYY+IY 
30 U(ND,2) = U (NDI, 2) 

PI = AT AN (1.0) *4.0 
DO 40 IX = NX1+2 , NX1+NX2 
DO 40 IY = 2, NY 
ND = (IX-1) *NYY*NZZ+NYY+IY 
NDI= NYY+IY 

U(ND,2) = U (NDI, 2) *COS ( (IX-NX1-1) *PI/NX2) 

40 U(ND, 3) =-U (NDI, 2) *SIN ( (IX-NX1-1) *PI/NX2) 

DO 45 IX = NX1+NX2+1,NX+1 

DO 45 IY = 2, NY 

ND = (IX-1) *NYY*NZZ+NYY+IY 

NDI= NYY+IY 

U (ND, 2) =-U (NDI, 2) 

45 CONTINUE 
C 

C U(ND, 3) AND U(ND, 4) ARE ZERO (HENCE, U(ND,5) IS AS DEFINED BELOW) 



oooooo 


DO 50 ND=1,NNM 
U (ND, 1 ) =DNST0 
U (ND,2) =U (ND, 2) *U(ND,1) 

U (ND, 6) =U (ND, 1) *R0*TEMP0 

U(ND, 5) =U (ND, 6) /GAM1+0 . 5*U (ND, 2) *U (ND,2) /U(ND, 1) 
50 CONTINUE 
RETURN 
END 


SUBROUTINE INVDET (A, B, DET) 

IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION A(3, 3) ,B (3, 3) 

G(Z1, Z2, Z3, Z4) =Z1*Z2-Z3*Z4 
F(Z1,Z2,Z3,Z4)=G(Z1,Z2,Z3,Z4)/DET 

C1=G (A (2, 2) ,A(3,3) ,A(2,3) ,A(3,2) ) 

C2=G (A (2, 3) ,A(3, 1) ,A(2, 1) ,A(3, 3) ) 

C3=G (A(2, 1) ,A(3,2) , A (2, 2 ) , A (3, 1) ) 
DET=A (1, 1) *C1+A (1,2) *C2+A(1, 3) *C3 
B(1,1)=F (A (2, 2) , A (3, 3), A (3, 2), A (2, 3) ) 

B (1, 2) =-F (A (1,2) , A (3, 3) ,A(1,3) , A (3, 2) ) 
B(l, 3)=F(A(1,2) , A(2, 3) ,A(1,3) , A(2,2) ) 

B (2, 1) =-F (A (2, 1) , A (3, 3) ,A(2, 3) , A (3, 1) ) 
B (2, 2) =F (A (1, 1) , A(3, 3) ,A(3, 1) , A(l, 3) ) 

B (2, 3) =-F (A (1, 1) ,A(2,3) , A (1, 3) ,A(2, 1) ) 
B(3, 1)=F (A (2, 1) ,A(3,2) ,A(3, 1) ,A(2,2) ) 

B (3, 2) =-F (A ( 1, 1) ,A(3,2) , A(l,2) ,A(3,1)) 
B (3, 3) =F (A (1, 1) ,A(2,2) ,A(2, 1) ,A(1,2) ) 
RETURN 
END 


SUBROUTINE MATMUL (A, B, C, M, N, L) 

IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION A (M, N) , B (N, L) , C (M, L) 

DO 20 1=1, M 

DO 20 J=l, L 

SUM=0 . 0 

DO 10 K=l, N 

10 SUM=SUM+A ( I , K) *B (K, J) 

20 C (I, J) =SUM 
RETURN 
END 


SUBROUTINE SHAPEL (XYZ, SF, DF, NDIM, NPE) 


SHAPE FUNCTIONS FOR LINEAR, ISOPARAMETRIC 3-DIMENSIONAL ELEMENT 
THIS SUBROUTINE EVALUATES THE SHAPE FUNCTIONS AND THEIR FIRST 
DERIVATIVES AT THE GAUSSIAN POINT XYZ 


IMPLICIT REAL* 8 (A-H,0-Z) 

DIMENSION XNODE (8, 3) , XYZ (NDIM) , SF(NPE) , DF (NDIM, NPE) 

DATA XNODE/ -1 . 0D0, 2*1 . 0D0, 2*-l . 0D0, 2*1 . 0D0, -1 . 0D0, 2*-l . 0D0, 2*1 0D0 
1, 2*-l . 0D0, 2*1. 0D0, 4*-l. 0D0, 4*1 . 0D0/ 

FCK (A, B, C) =0 . 125*A*B*C 
DO 20 1=1, NPE 
XNP 1=XYZ (1) * XNODE ( 1 , 1)+1.0 
YNP1=XYZ (2 ) *XNODE (I, 2) +1 . 0 
ZNP1=XYZ (3) * XNODE ( I , 3) +1 . 0 
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SF (I) =FCK (XNP1, YNP1, ZNP1) 

DF (1, I) =FCK (XNODE (1,1), YNP1, ZNP1) 
DF ( 2 , I ) =FCK ( XNP 1 , XNODE (1,2) , ZNP 1 ) 
20 DF (3,1) =FCK ( XNP 1 , YNP 1 , XNODE (1,3) ) 
RETURN 
END 


SUBROUT INE SURFGM ( K 1 , KG 1 , ELXY2 , DSURF , GNORM, NBS , NGP , NPE , ND IM ) 

GNORM ( I , J, K, L) I-TH COMPONENT OF 'NORMAL*DS' ON J-TH BOUNDARY 

SURFACE AT (K, L) GAUSS POINT 


IMPLICIT REAL*8 (A-H,0-Z) 

DIMENSION ELXYZ (NPE, NDIM) , DSURF (NDIM, NPE, 6, NGP, NGP) ,GC(3,3) , 
* GNORM (NDIM, NBS , NGP , NGP ) 

K0= (Kl+1 ) 12 
K2=K0+1 

IF ( K2 . EQ . 4 ) K2=l 
K3=K2+1 

IF (K2 . EQ . 3) K3=l 
DO 200 NGP 1=1, NGP 
DO 200 NGPK=1, NGP 


CALL GCSURF (GC, DSURF, ELXYZ, NGP I, NGPK, NGP, Kl, NDIM, NPE) 
DO 100 1=1, NDIM 
11 = 1+1 

IF ( II . EQ . 4) 11=1 
12 = 11+1 

IF (II .EQ.3) 12=1 


100 GNORM ( I, KG1, NGP I, NGPK) = (GC (II, K2) *GC (12, K3) -GC(I1,K3) *GC(I2, K2) 
1 *(-l)**Kl 


200 CONTINUE 
RETURN 
END 






SUBROUTINE TADMSH (X, Y, Z , IBNDC, KELSUR, NOD, NSURF, NNM, NBS, NDF, NEM 
* NPE) 

* 

* 

* 


★ 

★ 

* 

* 

* 

★ 

* 

★ 

★ 

★ 

* 

★ 

★ 

* 

★ 

★ 

★ 

* 

★ 
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PURPOSE 


MESH GENERATOR FOR TURN AROUND DUCT. 

TO GENERATE A THREE DIMENSIONAL MESH FOR A TURN AROUND 
DUCT. THE ELEMENT LIBRARY HAS THREE TYPES OF ELEMENTS 
VIZ. 8-NODED, 20 NODED, AND 27 NODED BRICK ELEMENTS. 
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* 1 / 1 / 

* 6 0 0 7 

★ 


* LINEAR RECTANGULAR ELEMENT 

★ 

* LIST OF VARIABLES : 

★ 

* 


* 

NX1 

= 

NUMBER 

OF 

DIVISIONS 

IN 

FLOW DIRN. IN 

PART 

1 (INLET) 

* 

NX2 


NUMBER 

OF 

DIVISIONS 

IN 

FLOW DIRN. IN 

PART 

2 (CURVE) 

★ 

NX 3 

- 

NUMBER 

OF 

DIVISIONS 

IN 

FLOW DIRN. IN 

PART 

3 (OUTLET) 

★ 

NY 

* 

NUMBER 

OF 

DIVISIONS 

IN 

RADIAL DIRECTION: 


★ 

NZ 

= 

NUMBER 

OF 

DIVISIONS 

IN 

Z-DIRECTION: 



★ 

NPE 

* 

NODES PER 

ELEMENT. (8 

OR. 20 OR 27) 




NOD (NNM, NPE) = CONNECTIVITY MATRIX 

- ELEMENT TYPE (1 - LINEAR (8 NODED) ; 2 = QUADRATIC) 
= INNER RADIUS OF THE CURVE. 

= OUTER RADIUS OF THE CURVE. 

= X - COORD INARE OF FIRST NODE IN X-Y-Z PLANE. 

= Y - COORDINARE OF FIRST NODE IN X-Y-Z PLANE. 

= Z - COORDINATE OF FIRST NODE IN X-Y-Z PLANE 

= ARRAY CONTAINING X-COORDINATES OF NODES. 

= ARRAY CONTAINING Y-COORDINATES OF NODES. 

= ARRAY CONTAINING Z-COORDINATES OF NODES. 


IEL 

R1 

R2 

XO 

YO 

ZO 

X (NNM) 
Y (NNM) 
Z (NNM) 




IMPLICIT REAL* 8 (A-H,0-Z) 

COMMON/MSH /ARCANG, NX, NY, NZ , NX1 , NX2 , NX3 

DIMENSION DX1 (10) ,DX3(10) ,DY(20) , DZ ( 5) , X (NNM) , Y (NNM) , Z (NNM) , 

# IBNDC (NNM, NDF) , KELSUR (NBS, 2) , NOD (NEM,NPE) 

READ (5,*) NX1, NX2, NX3, NY, NZ, IEL, NPE, Rl, R2, XO, YO, ZO, 

# ARCANG 

COMPUTE THE NUMBER OF ELEMENTS AND NODES IN THE MESH: 

PI = 3.141592654 

NELM = (NX1+NX2 +NX3) *NY*NZ 

NXX1 = IEL*NX1 

NXX3 = IEL*NX3 

NYY = IEL*NY 

NZZ = IEL*NZ 

IF (ZO .LE. 1.0E-10) THEN 
PHI = 0.0D0 

ELSE 

PHI = ATAN (Z0/X0) 

END IF 

ARCANG = ARCANG*PI/180 

ANGINC = ARCANG/NZZ 

RZ = DSQRT (Y0*Y0 + Z0*Z0) 

RZ = YO 

READ (5,*) (DX1 ( I ) , 1=1, NXX1 ) 

READ (5, *) (DX3 (I) , 1=1, NXX3) 

READ (5, *) (DY (I) , 1=1, NYY) 

NXX1 = IEL*NX1 + 1 
NXX2 = IEL*NX2 
NXX3 = IEL*NX3 
NYY = IEL*NY + 1 
NZZ = IEL*NZ + 1 


IF (NPE .EQ. 20) THEN 



NDS = NYY* ( (NX1 + NX2 + NX3 + 1)*(NZ+1)) + 

# (NY+1) * ( (NX1+NX2+NX3+1) + (NZ+1) * (NX1+NX2+NX3) ) 

ELSE 

NDS = NYY* (NXX1 + NXX2 + NXX3) *NZZ 
END IF 

IF (NDS .NE.NNM .OR. NEM . NE . NELM) THEN 

WRITE (6, 99 9 )NNM, NDS, NEM, NELM 

STOP 

ENDIF 

NTX = IEL*NX1 + 1 
NT XX = IEL*NX2 
NT XXX = IEL*NX3 
NTXT = NTX + NT XX 
NTXTT = NTXT + NTXXX 

COMPUTE THE NODAL COORDINATES IN SECTION 1 (STRAIGHT INLET) 

NTY = IEL*NY + 1 
NTZ = IEL*NZ + 1 
NY1 = ( IEL-1 ) *NY + 1 


IIX = 0 
L = 0 

DO 1050 IX = 1, NTX 

IF (NPE .EQ. 20) THEN 
MODY = MOD (IX, 2) 

ELSE 

MODY = 1 
END IF 


ZC = Z0 
ANGLE = PHI 

IF (MODY .EQ. 1) THEN 

IF (NPE .EQ. 20) THEN 

I = (NYY* (NZ+1) + (NY+1) * (NZZ) ) *IIX 

ELSE 

I = NYY* (IX - 1) *NZZ 
END IF 

DO 1020 IZ = 1, NTZ 

IF (NPE .EQ. 20) THEN 
MODZ = MOD ( IZ, 2) 

ELSE 

MODZ = 1 
END IF 

IF (MODZ .EQ. 1) THEN 
1 = 1 + 1 
X(I) = X0 

Y(I) = RZ*COS (ANGLE) 

2(1) = RZ*SIN (ANGLE) 

DO 1000 IY = 1, NTY-1 
1 = 1 + 1 
X(I) = X0 

Y ( I ) = (Y(I-l) + DY(IY) ) *COS (ANGLE) 
Z(I) = (Y(I-l) + DY(IY) ) *SIN (ANGLE) 
CONTINUE 


ELSE 

1 = 1 + 1 
X(I) = X0 

Y ( I ) = Y0*COS (ANGLE) 



Z(I) = ZC* S IN (ANGLE ) 


★ 

DO 1010 IY = 1, (NTY-NY1) 

1 = 1 + 1 
K - 2*IY - 1 
X(I) = X0 

Y (I) = (RZ + DY (K) + DY(K+1) ) *COS (ANGLE) 
Z<I) = (RZ + DY (K) + DY (K+l) ) *SIN (ANGLE) 
CONTINUE 

END IF 

IF ( IZ .LT. NTZ) ANGLE = ANGLE + ANGINC 
1020 CONTINUE 

IIX = IIX + 1 

ELSE 

DO 1040 IZ = 1, (NZ+1 ) 

1 = 1 + 1 
M = 2*IZ - 1 
X(I) = X0 

Y (I) = RZ*COS (ANGLE) 

Z(I) = RZ*SIN (ANGLE) 

DO 1030 IY = 1, (NTY-NY1) 

1 = 1 + 1 
K = 2*IX - 1 
X(I) = X0 

Y ( I ) = (RZ + DY (K) + DY (K+l) ) *COS (ANGLE) 

Z(I) = (RZ + DY (K) + DY (K+l) )*SIN (ANGLE) 
CONTINUE 

ANGLE = ANGLE + ANGINC 
CONTINUE 

END IF 

IF (IX .LT. NTX) X0 = X0 - DX1 (IX) 

1050 CONTINUE 
★ 

* COMPUTE THE NODAL COORDINATES IN THE CURVED SECTION" 

* 

NXPT1 = NTX + 1 
THINC = PI/NXX2 
THETA « PI + THINC 
YC = Y0 + R2 

* 

DO 1110 IX = NXPT1 r NT XT 
IF {NPE .EQ. 20) THEN 
MODY = MOD (IX, 2) 

ELSE 

MODY = 1 
END IF 

★ 

ZC = zo 

ANGLE = PHI 

IF (MODY .EQ. 1) THEN 

DO 1080 IZ = 1, NTZ 

IF (NPE .EQ. 20) THEN 
MODZ - MOD ( I Z , 2) 

ELSE 

MODZ - 1 
END IF 

* 

IF (MODZ .EQ. 1) THEN 
1 = 1 + 1 


1030 


1040 

★ 


1010 

★ 



* 


1060 

★ 

1070 

1080 

★ 

* 


★ 


1090 

1100 

* 

* 

1110 

* 

* 

* 

* 

* 


X(I) = X0 + R2*SIN (THETA) 

Y(I) = (YC + R2*C0S (THETA) ) *COS( ANGLE) 

Z(I) = (YC + R2*COS (THETA) )* SIN (ANGLE) 

DYY = 0.0D0 
DO 1060 IY = 1, NTY-1 
1 = 1 + 1 

DYY = DYY + DY (IY) 

X(I) = X0 + (R2 - DYY) *SIN (THETA) 

Y (I) = (YC+ (R2-DYY) *COS (THETA) ) *COS (ANGLE) 
Z (I) = (YC+ (R2-DYY) *COS (THETA) ) *SIN (ANGLE) 
CONTINUE 

ELSE 

1 = 1 + 1 

X(I) = X0 + R2*SIN (THETA) 

Y(I) = (YC + R2*COS (THETA) ) *COS (ANGLE) 

Z(I) = (YC + R2*COS (THETA) )* SIN (ANGLE) 

DYY = 0.0D0 

DO 1070 IY = 1, (NTY-NY1) 

1 = 1 + 1 
K = 2*IY - 1 

DYY = DYY + DY (K) + DY(K+1) 

X(I) = X0 + (R2 - DYY) *S IN (THETA) 

Y (I) = (YC+(R2-DYY) *COS (THETA) ) *COS (ANGLE) 
Z(D “ (YC+ (R2-DYY) *COS (THETA) ) *COS (ANGLE) 
CONTINUE 
END IF 

IF(IZ .LT. NTZ) ANGLE = ANGLE + ANGINC 
CONTINUE 
IIX = IIX + 1 


ELSE 


DO 1100 IZ = 1, (NZ+1 ) 

1 = 1 + 1 
M = 2*IZ - 1 

X (I) = X0 + R2*SIN (THETA) 

Y (I ) “ (YC + R2*COS (THETA) ) *COS (ANGLE) 
Z(I) = (YC + R2*COS (THETA) ) *SIN (ANGLE) 


DYY = 0.0D0 

DO 1090 IY = 1, (NTY-NY1) 
1=1 + 1 


K = 2*IY - 1 

DYY = DYY + DY (K) + DY(K+1) 

X(I) = X0 + (R2 - DYY) *SIN (THETA) 

Y(I) = (YC+ (R2-DYY) *COS (THETA) ) *COS (ANGLE) 

Z (I) = (YC+ (R2-DYY) *COS (THETA) ) *SIN (ANGLE) 
CONTINUE ' 


ANGLE = ANGLE + 2.0*ANGINC 
CONTINUE 


END IF 

THETA = THETA + THINC 
CONTINUE 


COMPUTE THE NODAL COORDINATES IN SECTION 3 (STRAIGHT OUTLET) 

NTXP11 = NTXT + 1 

Y0 = Y0 + 2 . 0*R2*COS (PHI) 

J = 0 


DO 1170 IX = NTXP11, NTXTT 



IF (NPE .EQ. 20) THEN 
MODY = MOD (IX, 2) 

ELSE 

MODY = 1 
END IF 

* 

J = J + 1 

X0 = XO + DX3 ( J) 

ZC = ZO + 2 . 0*R2*SIN (PHI) 

ANGLE * PHI 

★ 

IF (MODY .EQ. 1) THEN 

DO 1140 IZ = 1, NTZ 

* 

IF (NPE .EQ. 20) THEN 
MODZ = MOD (IZ, 2) 

ELSE 

MODZ = 1 
END IF 

* 

IF (MODZ .EQ. 1) THEN 
1 = 1 + 1 
X(I) = XO 

Y(I) = YO*COS (ANGLE) 

Z(I) = Y0*SIN (ANGLE) 

★ 

DYY = 0.0D0 

DO 1120 IY = 1, NTY-1 

DYY = DYY + DY (IY) 

1 = 1 + 1 
X(I) = XO 

Y (I) = (RZ + 2*R 2 -DYY) *COS (ANGLE) 

Z(I) = (RZ + 2*R2 - DYY) *SIN (ANGLE) 

1120 CONTINUE 

★ 

ELSE 

1 = 1 + 1 
X(I) = X0 

Y (I) = Y0*COS (ANGLE) 

Z(I) = Y0*SIN (ANGLE) 

★ 

DO 1130 IY = 1, (NTY-NY1 ) 

1 = 1 + 1 
K = 2*IY - 1 
X(I) = X0 

Y (I) = (RZ+2*R2-DY (K) -DY (K+l) ) *COS (ANGLE) 
Z(I) = (RZ+2*R2-DY (K) -DY (K+l) ) *SIN (ANGLE) 
1130 CONTINUE 

* 

END IF 

Hr 

IF ( IZ .LT. NTZ) ANGLE = ANGLE + ANGINC 
1140 CONTINUE 

IIX = IIX + 1 


ELSE 

* 


DO 1160 IZ = 1, (NZ+1 ) 

1 = 1 + 1 
M = 2*IZ - 1 
X(I) = X0*COS (ANGLE) 
Y (I) = Y0 

Z (I) = X0*SIN (ANGLE) 


DO 1150 IY = 1, (NTY-NY1) 
1 = 1 + 1 
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K = 2*IY - 1 
X(I) = XO 

Y(I) = (RZ+2*R2-DY (K) -DY (K+l) ) *COS (ANGLE) 
Z(I) = (RZ+2*R2-DY (K) -DY (K+l) ) *SIN (ANGLE) 
1150 CONTINUE 

ANGLE = ANGLE + 2 . *ANGINC 
1160 CONTINUE 

* 

END IF 

* 

1170 CONTINUE 
* 

DO 1175 1=1, NNM 
X(I)=0.0254*X(I) 

Y (I) =0 . 0254 *Y (I) 

1175 Z(I)=-0.0254*Z(I) 

* 

* DETERMINE THE CONNECTIVITY MATRIX: 

* 

NX = NX1 + NX2 + NX3 

IF (NPE .EQ. 20) NTY = 3*NY + 2 

★ 

DO 1200 IX = 1, NX 

DO 1190 IZ = 1, NZ 

DO 1180 IY = 1, NY 

* 

I = IY + (IX-1) *NY*NZ + (IZ-1) *NY 


IF (NPE .EQ. 20) THEN 

NOD (1,1) = IEL*IY - (IEL-1) + (NYY* (NZ+1) + 

(NY+1 ) *NZZ) * (IX-1) + (IZ-1) * (NYY+NY) 
NOD (1,2) = NYY* (NZ+1) + (NY+1) *NZZ + NOD (1,1) 
ELSE ' ‘ 

NOD (1,1) = IEL*IY - (IEL-1) + (IX-1)* 

(NYY*NZZ) *IEL+ (IZ-1) *IEL*NYY 
NOD (I, 2) = NYY*NZZ*IEL + NOD (1,1) 

END IF 

NOD (I, 3) = NOD (I, 2) + IEL 
NOD (I, 4) = NOD (1,1) + IEL 


* 


# 


# 


IF (NPE .EQ. 20) THEN 

NOD (1,5) = NTY + NOD (1,1) 

NOD (I, 6) = NYY* (NZ+1) + (NY+1) *NZZ + NOD (I, 5) 

ELSE 

NOD (I, 5) = NYY + NOD (1,1) 

NOD (I, 6) = NYY*NZZ*IEL + NOD (I, 5) 

END IF 


NOD (I, 7) « NOD (I, 6) + IEL 
NOD (I, 8) - NOD (I, 5) + IEL 
IF (NPE .EQ. 20) THEN 

NOD (I, 9) = NOD (1,1) + 
+ (1-IY) 


NYY* (NZ+1) + (NY+1 ) *NZ 


NOD (I, 10) 
NOD (I, 11) 
NOD (I, 12) 
NOD (I, 13) 
NOD (I, 14) 
NOD (I, 15) 
NOD (I, 16) 
NOD (I, 17) 

NOD (I, 18) 
NOD (I, 19) 
NOD (I, 20) 


NOD (1,2) + 1 
NOD (I, 9) + 1 
NOD (1,1) + 1 
NYY + NOD (1,1) 
NYY + NOD (I, 2) 
NOD (I, 14) + 1 
NOD (I, 13) + 1 
NOD (I, 5) + (NYY 
(1“IY) 

NOD (I, 6) + 1 
NOD (I, 17) + 1 
NOD (I, 5) + 1 


+ NY + 1) *NZ + 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

1180 

1190 

1200 

* 

* 

★ 


* 


ELSE IF (NPE .EQ. 
NOD (1,9) = 
NOD (1,10) = 
NOD ( I , 11) = 
NOD (I, 12) = 
NOD (I, 13) = 
NOD { I , 14) = 
NOD (1,15) - 
NOD (1,16) = 
NOD (I, 17) « 
NOD (I, 18) - 
NOD (I, 19) = 
NOD (1,20) - 
NOD (1,21) = 
NOD (1,22) - 
NOD (1,23) - 
NOD (1,24) - 
NOD (I, 25) - 
NOD (1,26) - 
NOD (1,27) * 
END IF 
CONTINUE 
CONTINUE 
CONTINUE 


27) THEN 

NOD (1,5) + NYY 

NOD (1,9) + NYY*NZZ*IEL 

NOD (I, 10) + IEL 

NOD (I, 9) + IEL 

NOD (1,1) + NYY*NZZ 

NOD (I, 2) + 1 

NOD (I, 13) + 2 

NOD (1,1) + 1 

NOD (I, 5) + NYY*NZZ 

NOD (I, 6) + 1 

NOD (1, 17) + 2 

NOD (I, 5) + 1 

NOD (I, 9) + NYY*NZZ 

NOD (I, 10) + 1 

NOD (I, 21) + 2 

NOD (I, 9) + 1 

NOD (I, 13) + 1 

NOD (I, 17) + 1 

NOD (I, 21) + 1 


COMPUTE THE NUMBER OF BOUNDARY SURFACES AND DETERMINE SURFACE 
INDICES 


NSURF = 2*NX* (NY+NZ) +2*NY*NZ 

* 

* ELEMENT FLUX SURFACES AT THE INLET OF THE DUCT: 

1 = 0 

NYZ = NY*NZ 
DO 1210 IYZ = 1, NYZ 
1 = 1 + 1 
KELSUR (1,1) = IYZ 
KELSUR (1,2) = 1 
1210 CONTINUE 

* 

* ELEMENT FLUX SURFACES AT THE SOLID SURFACE OF THE DUCT (OUTER) 

DO 1220 IX = 1, NX 

DO 1220 IZ = 1, NZ 
1 = 1+1 

ILL = (IX-1) *NY*NZ + (IZ-1) *NY + 1 
KELSUR (1,1) = ILL 
KELSUR (1,2) = 3 
1220 CONTINUE 

★ 

* ELEMENT FLUX SURFACES AT THE SOLID SURFACE OF THE DUCT (INNER) ; 

DO 1230 IX-1, NX 

DO 1230 IZ = 1, NZ 
1 = 1 + 1 

ILL = (IX-1) *NY*NZ + IZ*NY 
KELSUR (1,1) = ILL 
KELSUR (I, 2) = 4 
1230 CONTINUE 

* 

* ELEMENT FLUX SURFACES AT SYMMETRY SURFACE OF THE DUCT (IZ = 1) : 

DO 1240 IX-1, NX 

DO 1240 IY = 1, NY 
1 = 1 + 1 



ILL = IY + (IX-1) *NY*NZ 
KELSUR (1,1) = ILL 
KELSUR ( 1/ 2 ) = 5 
1240 CONTINUE 
* 

* ELEMENT FLUX SURFACES AT SYMMETRY SURFACE OF THE DUCT <IZ = NZ) 

* 

DO 1250 IX = 1, NX 

DO 1250 IY = 1, NY 
1 = 1 + 1 

ILL = IY + (IX-1) *NY*NZ + NY* (NZ-1) 

KELSUR ( 1/ 1) = ILL 
KELSUR (I, 2) = 6 
1250 CONTINUE 
★ 

* ELEMENT FLUX SURFACES AT THE INLET OF THE DUCT: 

* 

J = 0 

DO 1260 IZ = 1, NZ 

DO 1260 IY = 1, NY 
J = J + 1 
1 = 1 + 1 

ILL = (NX-1) *NY*NZ + J 
KELSUR (1,1) = ILL 
KELSUR (I, 2) = 2 
1260 CONTINUE 


DETERMINE THE BOUNDARY CONDITIONS: 

NBNDC = 0 
ND = 0 

NXX = NX + 1 
NYY = NY + 1 
NZZ = NZ + 1 


DO 1212 1=1, NDS 



DO 1212 

J = 1, 5 


IBNDC (I, 

J) = 1 

1212 

CONTINUE 


★ 



★ 

* 

SPECIFY THE 

INLET 


DO 1280 ID = 

1, NYY 


DO 1270 

JD = 1, NZZ 


ND 

= ND + 1 


NBNDC = NBNDC + 1 
IBNDC (ND, 2 ) = 0 
IBNDC (ND, 3) = 0 
IBNDC (ND, 4) = 0 
IBNDC (ND, 5) = 0 
1270 CONTINUE 

1280 CONTINUE 


BOUNDARY DEGREES OF FREEDOM 


* SPECIFY THE SOLID-WALL BOUNDARY DEGREES OF FREEDOM 

* 

DO 1300 KD = 1, NX 

ND1 = (NYY*NZZ)*KD + 1 
DO 1290 JZ = 1, NZZ 

ND = ND1 + (JZ-1) *NYY 
NBNDC = NBNDC + 1 
IBNDC (ND, 2) = 0 
IBNDC (ND, 3) = 0 
IBNDC (ND, 4) = 0 
CC IBNDC (ND, 5) = 0 



o o 



NBNDC = NBNDC 

+ 1 


IBNDC (ND+NY, 2 ) 

= 0 


IBNDC (ND+NY, 3) 

= 0 


IBNDC (ND+NY, 4) 

= 0 

* o 
o 

IBNDC (ND+NY, 5) 

= 0 

1290 

* 

CONTINUE 


1300 

* 

CONTINUE 


* 

SPECIFY THE EXIT 

BOUNDARY DEGREES OF FREEDOM 


* 


NBD1 = NYY*NZZ*NX 
DO 1320 1=1, NZZ 

NBD = NBD1 + (1-1) *NYY 
DO 1310 J = 1, NYY 
NBD = NBD + 1 
NBNDC = NBNDC + 1 
IBNDC (NBD, 5) = 0 
IBNDC (NBD, 3) = 0 
IBNDC (NBD, 4) = 0 
1310 CONTINUE 

1320 CONTINUE 
C 

RETURN 

999 FORMAT (/, 5X, ' ***** THE PARAMETERS NNM AND NEM SENT FROM THE MAIN 

# DO NOT COINCIDE WITH THOSE GENERATED IN TADMSH *****',/, 5x, ' ***** 

# THE PROGRAM IS TERMINATED *****',/ , 5X, ' NNM, NDS, NEM, NELM =',4I5) 
END 

FLOW IN A TURN-AROUND-DUCT (15X8X2 MESH) 

1 1 02 100 05. 1.0 0.8 

1.79E-03 293.0 110.0 287.0 1.402 0.72 0.1 1.205 

3 8 4 8 2 1 8 1.0 3.0 33.0 9.0 0.0 2.0 

20.0 8.0 2.0 

0.5 1.5 8.0 20.0 

0.1 0.2 0.3 0.4 0.4 0.3 


0.2 0.1 



